* 


AD-A254  349 

■llllllll 

PL-TR-92-2086 


ANALYSIS  AND  MODELING  OF  THE  CRRES  SATELLITE  DATA 


J.  N.  Bass 
C  A.  Hein 
C  E.  Jordan 
R.  P.  Vancour 


Radex,  Inc. 

Three  Preston  Court 
Bedford,  MA  01730 


March  20,  1992 


DTIC 

ELECT E 

JUL301992 

i 


D 


ScientiSc  Report  No.  10 


Approved  for  public  release;  distribution  unlimited 


PHILLIPS  LABORATORY 

AIR  FORCE  SYSTEMS  COMMAND 

HANSCOM  AIR  FORCE  BASE,  MASSACHUSETTS  01731-5000 


92  7  28  047 


92-20404 

I  miia  Mim  _ _  *  * 


"This  technical  report  has  been  reviewed  and  is  approved  for  pviblication" 


EDWARD  C.  ROBINSON 
Contract  Manager 
Data  Analysis  Division 


his  report  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  is 
eleasable  to  the  National  Technical  Information  Service  (NTIS) . 


ualified  requestors  may  obtain  additional  copies  from  the  Defense  Technical 
nformation  Center.  All  others  should  apply  to  the  National  Technical 
nformation  Service. 


f  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing 
ist,  or  if  the  addressee  is  no  longer  employed  by  your  organization,  please 
otify  GL/IMA,  Hanscom  AFB,  MA  01731.  This  will  assist  us  in  maintaining  a 
urrent  mailing  list. 


o  not  return  copies  of  this  report  unless  contractual  obligations  or 
otices  on  a  specific  document  requires  that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704  0188 


Public  reporting  burden  for  this  coHeciion  of  information  is  estimated  to  average  t  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  ano  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  asoea  of  this 
collection  of  information,  including  suggestions  tor  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  information  Operations  ano  Reports.  U IS  Jefferson 
Oavts  Highway.  Suite  1204.  Arlington.  VA  22202-4302.  and  to  the  Office  of  ^aanagement  and  Budget.  Paperwork  Reduction  Project  (0704-01 88).  Washington.  DC  20S03. 


1.  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  ANO  SUBTITLE 


2.  REPORT  DATE 
20  March  1992 


Analysis  and  Modeling  of  the  CRRES  Satellite  Data 


6.  AUTHOR(S) 

J.  N.  Bass,  C  A.  Hein, 

C  E.  Jordan,  R.  P.  Vancour 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS(ES) 

RADEX,In& 

Three  Preston  Court 
Bedford,  MA  01730 


3.  REPORT  TYPE  ANO  DATES  COVERED 
Scientific  Report  No.  10 


S.  FUNDING  NUMBERS 
PE  62101F 

PR 7659  TAOS  WUAB 
ContractF1962S«9e0068 


S.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


RXR-92031 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  ANO  AODRESS<ES) 
Phillips  Laboratory 
Hanscom  AFB,  MA  01731-5000 

Contract  Manager  Edward  C  Robinson/GPD 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


PL-TR-92-2086 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 
Approved  for  Public  Release 
Distribution  Unlimited 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 

The  radiation  belts  consist  of  energetic  charged  particles  trapped  by  the  Earth’s 
magnetic  field  which  consists  of  an  internal  Geld  and  an  external  Geld.  The  external  Geld  dominates  the  non-dipole 
portions  of  the  internal  Geld  at  and  above  geosynchronous  altitude.  External  Geld  modeling  usually  begins  with 
models  of  the  various  current  systems  from  which  the  Geld  is  derived.  CRRES  data  provides  for  a  generation  of  new 
models  of  the  radiation  belts.  A  Strawman  proposal  was  developed  for  the  CRRES  static  radiation  belt  model 
Product  Associated  Data  Base.  Several  groups  of  models  were  studied  in  preparation  for  working  with  the  CRRES 
data.  The  models  and  software  to  run  them  constitute  a  suite  of  routines  divided  into  Gve  categories.  A  Bj^^^ 
survey  data  base  was  constructed  to  study  the  effects  of  magnetic  latitude  and  local  time  at  The  behavior  of 

the  PROTEL  instrument  on  CRRES  was  modeled  in  its  hostile  environment  to  evaluate  its  performance.  Predictions 
of  the  developed  Contamination  Code  compared  very  well  with  the  results  from  the  Harvard  Cyclotron  experiments. 
The  CRRES  Dosimeters  were  modeled  for  an  isotropic,  a  mirror  plane  and  the  sin^a  pitch  angle  distributions. 


14.  SUBJECT  TERMS 

Radiation  belt  models,  Magnetic  Geld  models,  Magnetospheric  dynamics. 
Energetic  particles,  CRRES  satellite 


15.  NUMBER  OF  PAGES 

110 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UncIassiGed 


SECURITY  CLASSIFICATION 
OF  THIS  PAGE 
UncIassiGed 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

UncIassiGed 


20.  LIMITATION  OF  ABSTRACT 


Unlimited 


MSN  7540-01 -280-5500 


Standard  ^orm  ^98  2-89) 


TABLE  OF  CONTENTS 


Section  Page 

1.0  STATIC  RADIATION  BELT  MODELING  STUDIES  .  1 

1.1  Modeling  in  the  Near-Geosynchronous  Region  .  1 

1.1.1  General  Considerations .  1 

1.1. 1.1  Model  Overview .  2 

1.1. 1.2  Model  L  Values  and  Drift  Orbits  .  4 

1.1. 1.3  Third  Adiabatic  Invariant .  6 

1.1.2  Magnetic  Field  Model  Selection .  11 

1.1. 2.1  Strawman  Models .  11 

1. 1.2.2  Existing  Models .  18 

1. 1.2.2. 1  Constant  B  Equatorial  Contours .  18 

1.1.2.2.2  Polar  Cap  Equatorial  Projections  .  18 

1.1.3  Dosage  Calculations .  18 

1.2  Product-associated  Data  Base  Design .  24 

1.2.1  Strawman  Proposal .  24 

1.2.1. 1  Spatial  Binning  .  24 

1.2. 1.2  Sorting  Geographic  Bins  into  L  Bins .  25 

1.2. 1.3  Pitch  Angle  Binning .  25 

1.2. 1.4  Format  .  26 

1.2.2  L  vs  Pitch  Angle .  26 

1.3  Software .  28 

1.3.1  Software  Library .  28 

1.3.2  Other  Developments .  45 

1.3.2. 1  Conversion  of  Observed  Pitch  Angles  to  Equatorial .  45 

1.. 3.2. 1.1  Subroutine  BFLDSM  .  46 

1.3.2.1.2  CRRES  Coordinate  Transformation  Software 

Library .  49 

1. 3.2.2  L  vs  Pitch  Angle  for  Ephemeris  .  50 

1. 3.2.3  Adiabatic  Invariants  and  Dipole-Equivalent  Invariants .  50 

1. 3.2.3. 1  Equivalent  Parameters  .  52 

1.3.2.3.2  Software .  52 

1.4  References .  55 

2.0  DYNAMIC  MAGNETIC  FIELD  MODELS .  57 

2.1  Models  Published  Before  1989  .  57 

2.1.1  Model  Derivation  and  Description .  58 

2. 1.1.1  Olson-Pfitzer  Tilt  Dependent  Model  .  58 

2. 1.1. 1.1  Model  Derivation .  59 

2.1. 1.2  Olson-Pfitzer  Dynamic  Model .  59 

2. 1.1.3  Mead-Fairfield  Model .  60 


iii 


Table  of  Contents  (cont’d) 

Section  Page 

2. 1.1. 4  Tsyganenko-Usmanov  and  Tsyganenko  Models  .  60 

2.1. 1.4.1  Ring  Current  Model .  61 

2. 1.1. 4.2  Tail  Current  Models .  61 

2.1. 1.4.3  Magnetopause  Current  Models  .  61 

2.1. 1.4.4  Model  Derivation .  62 

2. 1.1.5  Stern  Parabolic  Magnetopause-Stretched  Magnetopause  Model  ...  62 

2.1. 1.5.1  Model  Derivation .  62 

2.1.2  Computation  Times  .  63 

2.1.3  L  Parameter  Comparisons .  63 

2.1.4  Models  vs  Data  .  66 

2.1.4.1  Comparisons  of  Models  with  SCATHA  Data .  71 

2. 1.4.2  Midnight  Equatorial  Field  Strength  Depression  .  75 

2.1.5  Summary .  77 

2.2  Models  Published  in  1989  or  Later .  78 

2.2.1  Tsyganenko  1989  Model  .  78 

2.2.2  Hilmer- Voigt  Model .  78 

2.3  References .  81 

3.0  EPHEMERIS  STUDIES .  84 

3.1  Magnetic  Latitude  and  Local  Time  at  Bmin .  84 

3.2  B/Bq  vs  L  and  Magnetic  Latitude .  87 

4.0  INSTRUMENT  MODELING .  90 

4.1  PROTEL  Contamination  Code .  90 

4.1.1  Motivation .  90 

4.1.2  Methodology  .  91 

4. 1.2.1  Physical  Assumptions .  91 

4. 1.2.2  Janni  Range-Energy  Tables .  91 

4.1.2.3  Ray  Tracing  of  PROTEL  HEH  .  91 

4. 1.2.4  Monte  Carlo  Technique  .  92 

4.1.3  Results .  92 

4. 1.3.1  Harvard  Accelerator  Experiment .  92 

4. 1.3.2  Contamination  Counts .  .  .  95 

4.2  Dosimeter  .  ...  95 

4.2.1  Motivation .  95 

4.2.2  Methodology  .  96 

4.2.2. 1  Physical  Assumptions .  97 

4.2.2.2  Path  Length  Distribution  Computations  .  98 

4.2.3  Results .  100 

4.3  References .  102 


iv 


List  of  Figures 


Figure 


Pagg 


1.  Baseline  magnetosphere  model  overview . 

2.  Differences  between  equatorial  L-shell  value  computed  for  various  models  and  the  dipole  L- 

shell  value  Lj . 

3.  Equatorial  traces  of  drift  shells  of  particles  crossing  the  midnight  equator  of  6.6  Rg,  with 

pitch  angles  40®  and  90® . 

4.  Third  invariant  . 

5.  Measured  SCATHA  proton  distributions . 

6.  Comparison  of  measured  and  Mead-ITJ  model  magnetic  field  magnitudes  for  SCATHA 

satellite,  20  April  1979  . 

7.  Schulz-McNab  source-surface  model . 

8.  Same  as  Figure  6,  except  for  Schulz-TU  model  instead  of  Mead-TU  model . 

9.  Field  line  configuration  for  Mead-TU  model  in  noon-midnight  meridian  . 

10.  Same  as  Figure  9,  except  for  Schulz-TU  model . 

11.  Equatorial  constant  B  contours  for  40®  pitch  angle  particles  as  determined  from  the 

Hilmer- Voigt  model  . 

12.  Flux  predicted  for  the  CRRES  HEEF  2.0  to  2.5  MeV  electron  channel  based  on  the  NASA 

model  AE8MAX,  along  with  the  measured  data,  before  the  March  1991  storm . 

13.  Flux  predicted  fo'  the  CRRES  HEEF  2.0  to  2.5  MeV  electron  channel  based  on  the  NASA 

model  AE8MAX,  along  with  the  measured  data,  after  the  March  storm . 

14.  This  shows  the  dose  calculated  at  the  center  of  the  LoLet  Dome  1  electrons  summed  over 

10  orbits  at  a  time . 

15.  Deviation  of  L  along  internal  field  model  lines  of  force  . 

16.  Deviation  of  L  along  a  magnetic  field  line  of  force  crossing  the  midnight  equator  at 

1-5  Re  . 

17.  Same  as  Figure  16,  except  for  2  Rg . 

18.  Same  as  Figure  16,  except  for  3  Rg . 

19.  Same  as  Figure  16,  except  for  4  Rg . 

20.  Same  as  Figure  16,  except  for  5  Rg . 

21.  Same  as  Figure  16,  except  for  6  Rg . 

22.  Same  as  Figure  16,  except  for  7  Rg . 

23.  Same  as  Figure  16,  except  at  noon  instead  of  midnight . 

24.  Same  as  Figure  17,  except  at  noon  instead  of  midnight . 

25.  Same  as  Figure  18,  except  at  noon  instead  of  midnight . 

26.  Same  as  Figure  19,  except  at  noon  instead  of  midnight . 

27.  Same  as  Figure  20,  except  at  noon  instead  of  midnight . 

28.  Same  as  Figure  21,  except  at  noon  instead  of  midnight . 

29.  Same  as  Figure  22,  except  at  noon  instead  of  midnight . 

30.  Predicted  CRRES  orbital  parameters  . 

31.  Magnetic  field  line  passing  through  a  given  vector  location  r  and  crossing  the  Equator  at 


32.  Overview  of  software  for  computation  of  adiabatic  invariants  and  dipole  equivalent 
invariants . 


3 

5 

7 

8 
10 

13 

14 

15 

16 
17 

19 

20 

21 

23 

27 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

47 

53 


V 


List  of  Figures  (c»nt’d) 


Figure  Page 

33.  Contours  of  constant  AL  =  L-L^:  internal  magnetic  field;  =  3.2.  3.7,  4.2,  and  4.7 .  64 

34.  Same  as  Figure  33,  but  for  Lj  =  5.2,  5.7,  6.2,  and  6.7 .  65 

35.  Same  as  Figure  33,  but  for  internal  model  +  Olson-Pfrtzer  tilt  dependent  external  model  ...  67 

36.  Same  as  Figure  35,  but  for  Lj  =  5.2,  5.7, 6.2,  and  6.7 .  68 

37.  Noon  and  midnight  equatorial  AL  profiles  for  1C  =  0 .  69 

38.  Same  as  Figure  37,  but  for  20*  N  magnetic  latitude .  70 

39.  Measured  -  model  differences:  20  April  1979  .  72 

40.  Noon  magnetopause  standoff  distance:  20  April  1979  .  73 

41.  Same  as  Figure  39,  but  for  21  April  1979  .  74 

42.  Averaged  measured  and  model  midnight  equatorial  magnetic  field  depression  .  76 

43.  Sample  of  a  preliminary  comparison  of  these  models  to  CRRES  data .  79 

44.  The  angle  between  the  CRRES  magnetometer  measured  field  vectors  and  the  model  field 

vectors  for  part  of  26  August  1990  .  80 

45.  Solar  magnetic  (SM)  local  time  and  orbit  number  for  all  cases  of  BMIN  distance  greater 

than  7  Rg,  orbits  1  -  450  .  85 

46.  SM  latitude  and  radial  distance  for  all  cases  of  BMIN  distance  greater  than  6  Rc,  orbits 

1  .450  .  . . .  86 

47.  CRRES  B/Bq  vs.  Dwell  Chart  -  1  Year  Statistics  .  88 

48.  CRRES  B/Bg  vs.  Solar  Magnetic  Latitude  at  L  =  2.5  (Orbits  1-510) .  89 

49.  PROTEL  HEH  Response  Function  for  the  isotropic  case  for  channel  10 .  93 

50.  Response  Function  for  Mirror  Plane  Case,  for  tw'o  orientations  of  the  mirror  plane  with 

respect  to  the  PROTEL  electron  deflection  magnet .  94 

51.  Etetector  Geometry .  98 

52.  Path  Length  Distribution  for  Detector  1  for  A,  =  60  degrees  compared  with  Infinite  Slab 

Path  Length  Distribution  (smooth  curve) .  100 

53.  Path  Length  Distribution  for  Detector  4,  X  =  60  degrees .  101 


List  of  Tables 


Table  Page 

1.  Sample  of  the  model  dose  tables  produced  for  the  CRRES  SpaceRad  dosimeter .  22 

2.  External  Magnetic  Field  Models .  57 

3.  Vax  8650  Computation  l  imes  (sec/10000  calls)  .  63 

4.  Average  Values  for  the  Olson-Pfitzer  Dynamic  Model .  75 

5.  Fraction  of  PROTEL  HEH  Counts  Due  to  Contamination .  95 

6.  Physical  and  Geometric  Parameters  of  the  CRRES  Dosimeter  Detectors  .  97 

7.  Notation  .  97 


ACKNOWLEDGEMENTS 


The  work  described  in  this  report  required  the  involvement  and  guidance 
of  a  number  of  individuals  at  PL.  and  their  interest  and  encouragement 
is  gratefully  acknowledged. 

Sue  Gussenhoven  and  Gary  Mullen  of  the  Space  Particle  Environment 
Branch  have  been  involved  throughout  the  projects,  and  provided 
essential  support  and  the  opportunity  for  many  challenging  studies. 

Ed  Robinson  of  PL7GPD  both  initiated  and  coordinated  the  activities  as 
Contract  Monitor.  We  also  deeply  appreciate  the  invaluable  help  and 
continuing  interest  of  Bob  Mclnerney,  Director,  Data  Analysis  Division. 

Howard  Singer  of  the  Plasma  Fields  and  Waves  Branch  provided 
magnetometer  data.  Robert  Hilmer  of  Rice  University,  Space  Physics 
and  Astronomy  Dept..  Karl  Pfitzer  of  McDonnell  Douglas  Cbrp..  Space 
Systems  Division  and  David  Stern  of  NASA  Grddard  Space  Flight  Center 
have  contributed  model  codes  and  revisions  as  necessary.  Bob  Filz  of  the 
Space  Particle  Environment  Branch  gave  useful  suggestions  and  assistance 
for  development  of  the  Protel  contamination  code  and  Dosimeter 
modeling  software.  We  have  also  benefitted  from  helpful  discussions  with 
many  other  unnamed  participants  in  the  CRRES  program. 


BTIC  QUALITY  liL 


vii 


iooessloD  For 

NTIS  GKAirl 
DTIC  T.-\R 


:ted 


tSUn  •\!inoviticfld 
JniJtif  Ic.Ttlori. 


Dlotrltutlon/ 

Av.allj^bl.lity  Cooes 

Diat 

1 

Avail 

Spec 

Mjd/or 

lal 

k' 

□  □ 


1.0  STATIC  RADIATION  BELT  MODELING  STUDIES 


1.1  Modeling  in  the  Near-Geosvnchronous  Region 


Ilie  radiation  belts  consists  of  energetic  particles  trapped  by  the  Earth’s  magnetic  field,  which  is  the  sum 
of  two  parts:  that  which  we  will  call  the  "internal"  magnetic  field  because  it  is  due  to  sources  (currents) 
in  the  Earth’s  interior,  and  that  which  we  will  call  the  "external"  magnetic  field  because  it  is  due  to 
sources  outside  the  Earth’s  surface.  Existing  radiation  belt  models  incorporate  only  the  internal  field, 
consisting  of  a  dominant  dipole  with  higher  order  multipoles.  However  the  significance  of  the  external 
field  increases  with  altitude.  As  we  shall  see,  the  external  field  dominates  the  non-dipole  portions  of  the 
internal  field  at  and  above  geosynchronous  altitude.  The  external  field  introduces  a  significant  noon- 
midnight  asymmetry  not  present  in  the  internal  field.  This  results  in  noon-midnight  asymmetries  in  the 
particle  populations,  notably  in  their  pitch  angle  distributions.  Local  time  dependent  pitch  angle 
distributions  have  been  reported  by  West  (1979]  and  Mullen  and  Gussenkoven  [1983],  In  addition,  the 
external  field  contains  a  dependence  on  the  time-varying  tilt  of  the  Earth’s  magnetic  dipole  toward  the 
Earth-Sun  line.  Iherefore,  even  in  magnetically  quiet  times  the  total  magnetic  field  is  not  truly  static. 

Models  based  on  the  internal  field  alone  generally  are  adequate  for  engineering  purposes,  such  as  the 
study  of  radiation  damage  on  a  satellite.  This  is  because  the  largest  radiation  dosage  is  at  low  altitudes 
where  the  internal  field  dominates.  However,  models  based  on  the  total  (internal  -t-  external)  field  would 
be  useful  for  scientific  understanding  of  magnetospheric  dynamics,  and  for  the  in-situ  detection  of 
substorms.  'Ihe  diffuse  auroral  zone  in  the  ionosphere  maps  magnetically  to  the  near-geosynchronous 
region  at  the  equator,  'fhe  two  regions  are  coupled  by  region  2  field  aligned  currents.  Vaayliunaa 
[1970]  derived  a  representation  of  the  field  aligned  currents  in  isotropic  plasma  in  terms  of  the  equatorial 
plasma  pressure  gradient  and  magnetic  field  geometry.  Heinemann  (1990)  has  generalized  this  to 
anisotropic  plasma.  Mullen  and  Guaaenhoven  ( 1983)  have  observed  that  20-400  kev  ions  contribute  a 
substantial  fraction  to  the  total  pressure  in  the  near-geosynchronous  region.  Baker,  et.  al.  [1978]  and 
Baaa,  et.  al.  [1989|  have  observed  variations  in  electron  pitch  angle  distributions  associated  with 
substorm  growth  phase  and  onset.  Model  quiet  pitch  angle  distributions  might,  therefore,  provide  a 
useful  baseline  against  which  to  compare  observed  pitch  angle  distributions  to  detect  and  predict,  in  real 
time,  the  occurrence  of  substorms. 


1.1.1  General  Considerations 

Following  an  overview  of  possible  model  construction  and  usage,  we  compare  the  L  values  computed 
from  (internal  +  external)  field  models  with  values  computed  from  an  internal  field  model  alone.  This 
will  illustrate  the  points  made  above  concerning  the  significance  of  the  external  field  above 
geosynchronous  altitudes,  the  noon-midnight  asymmetry  it  introduces,  and  its  tilt  dependence.  It  will  be 
seen  that  the  internal  field  can  be  truncated  to  the  quadrupole  terms.  Then  we  investigate  the 
dependence  of  the  third  adiabatic  invariant  on  mirror  point  L.  mirror  point  B,  and  the  dipole  tilt.  If  the 
dependence  on  dipole  tilt  is  negligible,  the  particle  model  may  be  accurately  defined  with  the  energy  and 
mirror  point  L  and  B  as  independent  variables.  Finally,  we  summarize  the  simplifications  w'e  conclude 
are  justified  in  this  modeling. 


1 


1.1. 1.1  Model  Overview 


Figure  1  shows  an  overview  of  possible  model  construction  and  use.  We  neglect  any  non-inductive 
electric  field  and  assume  that  the  time  variations  in  the  magnetic  field  are  sufficiently  slow  to  insure 
constancy  of  the  three  adiabatic  invariants: 


|i  =  — ? —  sin^a  =  — ? — 

2mB  2mB^ 

J  =  2pl 

♦  =  J  BdA 

Here  p  (magnetic  moment),  J,  and  ♦  are  the  first,  second,  and  third  adiabatic  invariants,  respectively. 
The  momentum  of  the  particle  is  p,  its  mass  is  m,  and  its  initial  pitch  angle  is  a.  In  the  bounce  motion 
of  the  particle’s  drift  guiding  center  along  a  field  line,  is  the  magnetic  field  magnitude  at  the  mirror 
points,  where  the  pitch  angle  is  90° .  From  constancy  of  the  first  adiabatic  invariant,  this  is  given  by: 

=  -4-  (2) 

sin^a 


I  is  the  integral  of  the  bounce  motion  between  the  mirror  points: 


aN 


-±d. 


(3) 


where  A  and  A'  are  the  mirror  points  and  ds  is  an  element  of  length  along  the  field  line.  'Hie  mirror 
point  Mcllwain  L  value  is  a  function  of  and  I,  as  specified  by  Mcllwain  [1961],  Its  dependence 
on  the  pitch  angle  follows  from  the  dependence  of  L  on  the  mirror  points  in  a  non-dipole  field.  In  the 
third  adiabatic  invariant  the  integral  is  over  any  surface  enclosed  by  the  drift  shell  of  field  lines 
containing  the  instantaneous  locus  of  mirror  points  of  common  1  and  (or  and  B^). 

In  a  purely  static  field,  with  no  electric  field,  the  energy  and.  therefore,  the  momentum  p  are  constants 
of  the  motion.  Therefore,  it  follows  that  the  mirror  point  magnetic  field  B^,  the  bounce  integral  1,  and 
the  mirror  point  L  value  are  also  constants.  Ihe  drift  shell,  defined  by  I  and  B^  (or  and  B^), 
is  fixed  in  time.  From  this,  it  follows  that  constancy  of  the  third  adiabatic  invariant  is  not  an 
independent  condition.  'Thus,  it  seems  convenient  to  select  ener©f,  and  B^  as  the  independent 
variables  of  the  model.  In  this  case,  the  third  box  in  our  figure  may  be  eliminated,  and  the  variables  p, 
J  and  9  replaced  by  the  variables  L^,  and  E  in  subsequent  boxes. 
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BASELINE  MAGNETOSPHERE  MODEL  OVERVIEW 


Figure  I.  Baseline  magnetosphere  model  overview 
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In  a  time  varying  field,  the  third  invariant  in  general  depends  on  time  as  well  as  on  the  mirror  point 
parameters  of  the  drift  orbit.  Therefore,  it  is  necessary  to  change  these  mirror  point  parameters  in  order 
to  preserve  the  constancy  of  the  third  invariant.  This,  in  turn,  requires  a  change  of  the  particle’s 
momentum  and  energy  to  preserve  the  first  two  adiabatic  invariants.  Physically,  the  particles  are 
responding  to  an  induction  electric  field  due  to  the  time  derivative  of  the  magnetic  flux  enclosed  by  the 
drift  shell.  However,  if  the  dependence  of  the  third  invariant  were  sufficiently  "small”,  these  effects  may 
be  neglected:  the  third  invariant  would  be  approximately  constant  for  fixed  mirror  point  parameters, 
allowing  these  and  the  energy  to  remain  constants.  Physically,  the  variation  of  the  particle’s  energy  due 
to  time  variations  in  the  magnetic  field  would  be  negligible. 

If  the  time  variation  is  periodic,  as  is  the  case  for  the  dipole  tilt,  we  may  evaluate  its  maximum  effect. 
This  may  then  be  compared  with  other  sources  of  error,  for  instance,  measurement  error  and  magnetic 
field  model-to-model  differences. 


1.1. 1.2  Model  L  Values  and  Drift  Orbits 

Figure  2  illustrates,  for  four  magnetic  field  models,  the  difference  A  L  between  the  L  value  computed  with 
the  model  and  that  computed  with  a  pure  dipole  model.  Riese  are  plotted  vs  Lj.  the  pure  dipole  L. 
Each  mode!  is  the  sum  of  the  IGRF  1985  internal  field  model  and  the  indicated  external  field  model, 
except  for  the  model  designated  "No  Ext.",  which  is  the  IGRF  1985  alone.  ITe  external  models  used  in 
the  other  three  cases  are  the  I'syganenko  1989  model  \Tsyganenko,  1989],  the  Hilmer-Voigt  model 
[Hilmer,  1989|,  and  the  quiet  Olson-Pfitzer  model  \Olson  and  Pfitzer.  1977|.  Ihe  results  are  given  at 
the  equator  for  "Kp  =  0".  Actuall”,  only  the  Isyganenko  model  contains  explicit  Kp  dependena*.  The 
Hilmer-Voigt  model  is  the  "super-quiet"  Ciise,  while  the  Olson-Pfitzer  model  is  known  to  be  accurate  for 
quiet  conditions.  The  IGRF  1985  model  itself  contains  neither  explicit  nor  implicit  dependence  on 
magnetic  activity.  The  dipole  reference  was  computed  assuming  that  the  dipole  is  at  the  center  of  the 
Earth.  In  this  case  is  the  distance  in  R^  from  the  center  of  the  Earth.  In  all  of  the  plots,  the  sun’s 
magnetic  longitude  is  0°. 

I  he  IGRF  1985  results  (solid  lines)  indicate  the  contribution  of  the  non-dipole  terms  of  the  internal  field 
to  1  .  rtiis  is  clearly  smaller  than  the  contributions  due  to  the  external  field  I'he  IGRF  1985  AL  is 
nearly  constant  over  the  range  shown,  indicating  that  the  contribution  is  dominated  by  the  quadrupole 
terms.  Ihis  follows  from  the  fact  that  the  contribution  due  to  the  nth  degree  internal  terms  should  be: 


-(„-!))  ^  ^2-n 


(4) 


The  proportionality  follows  from  the  fact  that  the  strength  of  the  nth  degree  contribution  to  the  magnetic 
field  varies  with  distance  as  r  "  ^  Thus,  contribution  of  any  multipole  term  relative  to  the  dipole  (n=l) 
should  be  proportional  to  Lj  If,  instead  of  the  center  of  the  Earth,  we  used  the  position  of  the 
offset  (eccentric)  dipole  [Chapman  and  Bartels.  1940J  as  the  coordinate  origin,  AL  for  the  internal  field 
would  be  less  than  0.01  R^;  for  the  entire  range  of  the  plots. 
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Figure  2.  Differences  between  equatori;il  L-shell  value  computed  fur  various  models  and  tlie  dipole  L-shell  value  L 


Figure  3  shows  a  comparison  of  the  equatorial  crossings  of  the  drift  shells  for  particles  with  90°  and  40° 
equatorial  pitch  angles  at  geosynchronous  altitude  (6.6  R^)  at  local  midnight.  We  see  that  the  particles 
which  separate  at  local  times  other  than  midnight  (drift  shell  splitting),  for  the  equatorial  crossing  of  the 
40° ,  is  nearly  circular. 


1.1. 1.3  ITiird  Adiabatic  Invariant 

If,  for  given  values  of  a  drift  orbit’s  mirror  point  parameters,  the  third  adiabatic  invariant  is  independent 
of  time,  constancy  of  the  mirror  point  parameters  follows  from  constancy  of  the  adiabatic  invariants. 
Although  the  time-varying  dipole  tilt  rules  out  strict  time  independence,  the  variations  introduced  may 
be  sufficiently  small  that  constant  mirror  point  parameters  may  still  be  assumed.  Therefore,  it  is  useful 
to  study  the  tilt-dependence  of  the  third  invariant. 

Figure  4  shows,  on  the  left  hand  panel,  the  L  dependence  of  the  third  invariant  for  a  particle  with  zero 
second  adiabatic  invariant  (90°  equatorial  pitch  angle)  for  the  dipole  field,  and  for  the  dipole  -H 
Tsyganenko  1989  magnetic  field  model  at  the  lowest  Kp  level.  In  the  case  here  of  zero  second  adiabatic 
invariant: 

Bn,  (1-0)  -  ~  (S) 


where  M  is  the  Earth’s  dipole  moment.  In  the  present  case,  M  was  derived  from  the  IGRF  85  internal 
field  model,  updated  to  1990.  The  right-hand  panel  of  the  figure  shows  dipole  tilt  dependence  of  the 
third  invariant  for  L  =  5.  6,  7  and  8  (for  zero  second  invariant)  as  a  percent  difference  from  the  value 
at  zero  tilt.  These  are  plotted  over  the  range  of  allowable  positive  tilt  values  (the  figure  for  negative  tilt 
would  be  the  mirror  image). 

The  left-hand  panel  shows  that  the  external  field  contributes  a  nearly  constant  negative  offset,  lliis  is 
due  to  the  nature  of  the  AB  contribution  of  the  external  field  to  the  total  field.  At  small  distances, 
inside  all  the  drift  shells  plotted  here,  AB  is  large  and  negative  at  all  local  times.  At  larger  distances, 
where  the  drift  orbits  are  located.  AB  is  much  smaller  in  magnitude  and  becomes  positive  on  the  dayside 
as  we  approach  the  magnetopause,  while  remaining  negative  in  the  nightside  tail.  The  large  negative  AB 
region  is  inside  the  orbits  plotted  here,  making  the  same  contribution  to  the  flux  enclosed  by  all  of  them. 
In  the  region  where  the  orbits  lie,  AB  is  smaller  in  magnitude  than  the  region  inside,  with  some  tendency 
toward  cancellation  of  the  day  and  night  portions.  Thus,  the  third  invariant  for  these  orbits  is  dominated 
by  the  contributions  from  their  common  inner  region.  The  right-hand  panel  shows  that  the  maximum 
variation  magnitude  due  to  tilt  is  -2.5  %,  at  the  maximum  tilt.  The  variation  increases  with  L  except 
for  a  reversal  between  L=7  and  L=8. 

To  gauge  the  impact  of  the  dependence  of  the  third  invariant  on  dipole  tilt,  we  look  at  two  aspects:  the 
dependence  of  a  particle’s  L  and  energy  on  tilt  that  is  required  to  preserve  its  adiabatic  invariants,  and 
the  resulting  tilt  dependence  of  fluxes  at  fixed  energy  and  L  From  the  theory  of  small  perturbations: 
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THIRD  INVARIANT  -  DIPOLE  +  TSYGANENKO  1983 


difference  from  zero  tilt  value  for  dipole  plus  Tsyganenko  1989  model. 
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Equation  (6)  describes  the  adjustment  in  L  that  a  particle  must  make  to  compensate  for  a  change  in  the 
third  adiabatic  invariant  at  fixed  L  The  second  approximate  equality  follows  from  the  approximation 
that  the  partial  is  -1,  as  for  a  dipole,  which  seems  justified  by  the  left-hand  panel  of  Figure  4.  Equation 
(7)  follows  from  the  conservation  of  the  first  adiabatic  invariant  (non-relativistic  particles),  and  the 
relationship  (Equation  (5))  between  B  and  L  for  zero  second  invariant. 

Given  the  maximum  percent  variation  of  the  third  invariant  with  tilt,  -2.5  %.  we  find  that  the  maximum 
L  adjustment  of  a  particle  is  -2.5  %,  and  the  maximum  energy  adjustment  is  -7.5%.  The  L  adjustment 
is  within  the  model-to-model  differences  indicated  in  Figure  2.  The  energy  adjustment  is  within  typical 
energy  widths  of  the  SCATHA  particle  detectors  ,-100  %  [Hanaer,  et.  al.,  1979],  but  exceeds  those  of 
some  of  the  CRRES  detectors  [Guaaenhoven,  et.  al.,  1985]. 

Using  the  Liouville  Theorem,  we  estimate  the  variation  in  the  particle  distribution  function  value  at  fixed 
E  and  L  as; 


^  =  -^LMnL  -  ^!!^MnE  (8) 

/  dL  dlnE 


To  get  some  feel  of  the  partials  of  the  distribution  function  f,  we  examined  1  MeV  protons  observed  by 
the  SCATHA  SC5  instrument  on  20  April  1979.  The  left-hand  panel  of  Figure  5  shows  the  observed 
(symbols)  L  dependence  of  the  1  MeV  particle  fluxes.  They  fit  well  the  expression  (solid  line) 


/=0.045€-2*^<‘'-*) 


(9) 


At  L=7,  the  energy  dependence  of  the  proton  population  (right-hand  panel.  Figure  5)  fits  well  the 
expression: 
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Figure  5.  Measured  SCATHA  proton  distributions.  Left  panel:  radial  distribution  at  1.1  MeV. 
Right  panel:  energy  distribution  at  L  =  7. 
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Plugging  these  results  into  Equation  (8),  we  find  the  maximum  variation  in  particle  fluxes  at  fixed  L  and 
E  due  to  dipole  tilt  variation  is  -11  %. 


1.1.2  Magnetic  Field  Model  Selection 


1. 1.2.1  Strawman  Models 

Modeling  of  the  near-geosynchronous  region  has  been  attempted  previously  employing  a  Mead-type 
magnetic  model  (defined  by  only  two  coefficients)  adjusted  to  best  fit  the  measured  magnetic  field  for 
the  period  in  question  [RADEX,  Ine.,  1989].  lliis  model  has  the  advantage  that  the  adiabatic  invariants 
can  be  computed  from  simple  analytic  expressions.  The  calculations  were  based  on  the  Pennington-Stern 
formulation,  applicable  to  a  field  expressed  as  the  gradient  of  a  scaler  potential,  which  itself  is  expressed 
as  an  expansion  in  internal  and  external  spheriail  harmonics.  The  results  are  accurate  through  first  order 
in  the  coefficients.  Since  the  chosen  model  is  independent  of  dipole  tilt  and  other  time-varying  dynamic 
parameters,  it  is  a  truly  static  model.  In  this  case,  the  third  adiabatic  invariant  is  strictly  a  function  of 
the  first  two;  therefore,  we  need  be  concerned  only  with  conservation  of  the  first  two  adiabatic  invariants. 

A  major  problem  with  this  form  of  model  is  that  it  fails  to  account  properly  for  the  decay  of  the  ring 
current  portion  of  the  field  with  distance,  since  the  azimuthally  symmetric  part  of  the  field  in  this  model 
is  spatially  constant.  Another  possible  problem  with  the  original  model  is  the  lack  of  dependence  on 
the  dipole  tilt  relative  to  the  Earth-Sun  line.  Although,  as  we  have  discussed  previously,  the  third 
adiabatic  invariant  depends  only  moderately  on  the  tilt  at  given  values  of  the  first  two  invariants,  it  is  still 
necessary  to  represent  the  tilt  variation  in  the  field  for  accurate  determination  of  the  first  two  invariants. 
A  third  problem  with  the  original  model  is  that  the  internal  part  of  the  field  includes  only  a  dipole.  We 
found  in  our  magnetic  field  model  studies,  that  the  higher  order  terms  produce  a  radially  independent 
contribution  to  L  which  is  in  excess  of  0.1  in  the  region  of  the  South  Atlantic  Anomaly.  The  Earth’s 
rotation  produces  a  small  time  variation  which  could  impact  the  third  invariant. 

A  model  proposed  here  will  hopefully  remove  these  difficulties  while  retaining  some  of  the  simplifications 
of  the  Mead  model.  The  first  modification  is  that  the  ring  current  would  be  represented  by  a 
Tsyganenko  form  which  is  simple,  analytic,  azimuthally  symmetric,  and  has  the  proper  behavior  at  large 
distances.  A  sum  of  two  such  terms  of  opposite  directions,  a  la  Hilmer,  could  also  be  considered.  The 
second  modification  is  that  the  external  spherical  harmonics,  which  would  continue  to  represent  the  tail 
and  magnetopause  portions  of  the  field  would  contain  tilt-dependent  terms.  These  would  be  in  the  form 
of  harmonics  whose  contributions  to  the  z  component  of  the  field  are  antisymmetric  with  respect  to  the 
dipole  equator.  The  third  modification  would  be  the  addition  of  a  quadrupole  portion  of  the  internal 
field.  An  analysis  of  the  Pennington  equations  indicates  that  these  terms  produce  the  above-mentioned 
radially  independent  contribution  to  L,  while  higher  multipoles  yield  contributions  which  vary  as  inverse 
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powers  of  the  distance.  By  placing  the  origin  at  the  offset  dipole  rather  than  at  the  Earth’s  center,  we 
can  eliminate  some,  but  not  all  of  this  effect. 

The  calculations  of  the  second  and  third  adiabatic  invariants  may  be  done  with  the  simplification,  used 
in  the  earlier  effort,  that  only  first  order  contributions  need  be  considered.  'Ihus,  contributions  from  each 
part  of  the  field  can  be  considered  separately.  The  spherical  harmonic  contributions  can  be  obtained 
from  the  Pennington-Stern  theory,  as  before.  This  leaves  the  contribution  of  the  ring  current  to  the 
adiabatic  invariants  to  be  determined.  Ihis  can  be  done  by  analytically  fitting  or  numerically  tabulating 
the  invariants  computed  for  a  dipole  +  ring  current  field  magnetic  field,  and  subtracting  the  dipole 
portion.  Both  procedures  are  simplified  by  the  independence  of  the  ring  current  model  on  local  time, 
longitude,  and  dipole  tilt.  Furthermore,  if,  consistent  with  the  inclusion  of  first  order  perturbations  only, 
we  assume  that  the  contributions  are  linearly  proportional  to  the  field  strength  parameter  of  the  ring 
current,  then  we  need  only  fit  or  tabulate  the  dependence  on  magnetic  latitude,  longitude,  and  the  ring 
current  length  scale  parameter. 

Figure  6  shows  a  model  constructed  from  a  Tsyganenko-Usmanov  ring  current  field,  an  internal  field 
through  the  n=2  terms,  and  the  tvio  leading  terms  of  the  Mead  model,  for  day  110  (1979).  The  Mead 
terms  have  been  computed  for  the  time-varying  standoff  distance  on  that  day.  The  Tsyganenko  term  is 
for  length  scale  parameter  equal  to  4  Rj;,  and  for  magnetic  field  strength  parameter  -65  nT  so  that  the 
net  contribution  of  the  Ring  current  and  Mead  magnetopause  terms  to  the  field  at  the  Earth’s  surface 
is  -40  nT  when  the  standoff  distance  is  10  Rj;.  Also  shown  are  data  (symbols)  from  the  SCATHA 
magnetometer.  There  is  an  almost  constant  difference  between  the  data  and  the  model,  with  the  model 
field  consistently  higher  than  the  data.  We  would  expect  a  tail  field  to  reduce  this  disagreement,  since 
it  would  produce  a  negative  contribution  in  the  inner  magnetosphere. 

Therefore,  the  source-surface  model  of  Schulz  and  McNab  [1986]  has  been  explored.  This  model 
contains  an  external  spherical  harmonic  representation  of  the  tail  and  magnetopause  fields  which  would 
be  valid  within  the  inner  magnetospheric  region  S  (Figure  7,  originally  Figure  1  of  that  paper).  The  two 
leading  terms  are  identical  in  form  to  the  two  Mead  model  terms  mentioned  above.  However,  the 
coefficient  of  the  symmetric  term  is  smaller,  and  the  coefficient  of  the  asymmetric  term  is  larger  than  the 
coefficients,  used  above,  which  represent  only  the  magnetopause  contribution.  Ihis  results  in  an  average 
decrease  in  the  inner  magnetospheric  field,  with  a  day-night  asymmetry,  which  is  exactly  what  one  expects 
from  a  tail  field.  Figure  8  shows  a  comparison  of  the  data  with  this  revised  model,  where  the  ring 
current  amplitude  is  readjusted  to  again  give  a  total  field  depression  of  -40  nT  at  the  Earth’s  surface  for 
10  Rf;  standoff  distance.  We  see  that  there  is  much  better  agreement  between  this  model  ("Schulz-TU”) 
and  the  data,  than  between  the  previous  model  ("Mead-ITJ")  and  the  data.  In  Figures  9  and  10,  the  field 
line  configurations  are  plotted  for  these  models,  for  10  R^  standoff  distance  and  including  the  n=2 
internal  terms.  Inside  10  R^,  the  Schulz-TU  model  looks  quite  reasonable.  The  divergence  of  the  field 
lines  outside  this  region  does  not  correctly  represent  that  model.  A  different  formulation,  which 
smoothly  joins  the  field  of  the  inner  region  S  at  the  boundary  between  the  two  regions,  is  prescribed  in 
the  tail  region  T.  as  shown  in  Figure  7.  Ihe  determination  of  the  coefficients  by  Schulz  and  McNab  was 
in  the  absence  of  a  ring  current  field.  Coefficients  resulting  from  inclusion  of  such  a  field  might  improve 
the  agreement  with  the  data. 

The  results  here  are  for  zero  dipole  tilt,  which  is  what  Schulz  and  McNab  considered.  Considerations 
by  Stern  [1985]  indicate  that  the  field,  in  more  general  cases,  may  be  represented  by  the  form: 
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Comparison  of  measured  and  Mead-TU  model  magnetic  field  magnitudes  for  SCATHA  satellite,  20  April  1979. 


Figure  7.  Schulz-McNab  source-surface  model  (Figure  1  of  Schulz  and  McNab  [1986]). 
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Figure  9.  Field  line  configuration  for  Mead-lU  model  in  noon-midnight  meridian. 
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where  t  is  the  tilt  angle,  Bq  is  the  model  obtained  by  Schulz  and  McNab  for  zero  tilt,  and  B^  is  the 
model  they  would  obtain  for  90°  tilt.  Ilius,  solutions  would  have  to  be  obtained  only  lor  these  two 
cases,  the  results  of  which  could  be  superposed  to  produce  the  model  for  arbitrary  tilt. 


1.1. 2.2  Existing  Models 

1.1. 2.2.1  Gjnstant  B  Equatorial  Contours 

In  an  effort  to  get  a  better  handle  on  how  to  label,  as  well  as,  model  the  trapped  particles,  we  needed 
to  develop  a  better  understanding  of  the  correlation  between  particles  at  high  latitude  and  low  altitude 
with  those  at  the  equator.  To  this  end,  we  made  three  sets  of  i-.aps  of  constant  equatorial  B  contours 
for  90°,  40°  and  10°,  respectively.  Several  different  external  field  models  were  used  with  the  IGRF  1980 
internal  field  model,  namely.  Tsyganenko  1989,  Hilmer- Voigt,  Olson-Pfitzer  Tilt-Dependent,  Mead 
Asymmetric  and  Mead  Symmetric.  In  addition,  IGRF  1980  alone  and  the  dipole  were  used  for  the 
SCATHA  Day  110  (  20  April  1979)  scenario.  Plots  were  produced  of  the  B  contours  for  each  of  the 
models  (Figure  J 1)  and  for  the  B  contours  for  a  given  model  with  the  dipole  B  a>ntours  superposed. 

1.1. 2.2.2  Polar  Gip  Equatorial  Projections 

In  addition  to  the  equatorial  constant  B  contours  discussed  above,  there  was  also  interest  in  mapping 
particles  which  mirror  at  840  km  and  100  km  back  to  the  equator  and  determining  their  constant  B 
a)ntours  as  well.  ITe  software  was  setup  to  do  this  (only  slight  modifications  to  the  software  were 
required  from  the  previous  work),  however,  no  plots  were  ever  generated.  Due  to  other  more  pressing 
pre-launch  work,  this  study  was  set  aside. 


1.1.3  Dosage  Calculations 

Before  the  CRRES  launch,  we  were  interested  in  simulating  the  data  and  creating  some  mock  data  files 
for  several  of  the  instruments  (Figures  12  and  13),  including  the  SpaceRad  dosimeter.  In  order  to  do 
this,  we  slightly  modified  the  Yates  and  Rothwell  dosage  code  presently  in  u.se  here  at  PL  lliis  code 
takes  a  given  orbit,  cx)mputes  the  fluence  from  the  NASA  trapped  r.idiation  nu)dels,  and  then,  ailculates 
the  dose  deposited  at  the  center  of  an  aluminum  sphere.  We  modified  the  default  thickness  and  energy 
parameters  to  just  include  the  specifications  related  to  the  dosimeter.  Some  prelaunch  files  were 
produced.  Following  the  launch,  we  began  to  produce  the  model  doses  per  orbit  based  on  orbital 
parameters  from  the  ephemeris  files.  In  addition,  we  kept  a  cumulative  total  of  each  of  the  HiLE  f  and 
LoLET  channels  (Table  1).  Additionally,  we  produced  plots  of  10-orbit-sums  of  the  dose  (Figure  14). 
fhese  tables  and  plots  are  regularly  updated  as  new  ephemeris  files  become  available. 

In  addition  to  this  regular  work,  two  special  cases  for  dosages  were  calculated.  Ihe  dose/day  over  one 
week  of  the  DMSP  F7  orbit  was  calculated  for  both  solar  minimum  and  maximum  cc>nditions.  Doses 
were  also  calculated  using  the  AP8MAX  model  for  the  HiLET  channels  of  the  SpaceRad  dosimeter. 
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SIvI  Coordinates  (R^-) 


_j  Q  Equatorial  Constant  B  Cantours  Hilmer  (40)  2^  April  1979  IGRF80 


X  SM  Coordinates  (Rj.) 


Fiuure  11.  This  shows  the  equatorial  constant  B  contours  for  40»  pitch  angle  particles  as  determined 
from  the  Hilmer- Voigt  model.  Ilie  jagged  edge  is  where  the  model  is  cut  off  due  to  the  magnetopause. 


2.0  to  2.5  MeV  Electrons 


AILitude  (RE) 


Figure  12.  This  figure  shows  the  flux  predicted  for  the  CRRE.S  HEEF  2.0  to  2.5  MeV  electron  channel 
based  on  the  NASA  model  AE8MAX,  along  with  the  me.isured  data,  before  the  March  1991  storm,  lliis 
is  representative  of  the  types  of  simulations  which  were  run  for  various  electron  and  proton  instruments. 
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Dose  for  Solid  Sphere  (Rods/10  Orbits)  ♦dO’ 

60  7.0  8.0  9.0  10.0 


NASA  Model  Dose  Summed  Over  10  Orbits 

SpoceRod  Dosimeter,  Dome  1,  LoLET  Electrons,  1.0-1125  KW 


Figure  14.  This  shows  the  dose  calculated  at  the  center  of  the  LoLet  Dome  1  electrons  summed  over 
10  orbits  at  a  time. 
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Tables  of  these  doses  were  based  on  one  second  fluences  and  were  binned  at  1/20  L-shell  fix)m  L*  1  to 
3  Rg  for  3  different  values  of  1.00, 1.35  and  2.24. 


1.2  Product-associated  Data  Base  Desian 

1.2.1  Strawman  Proposal 

A  strawman  proposal  was  developed  for  the  CRRES  static  radiation  belt  model  product-associated  data 
base.  This  concentrated  on  the  structure,  including  the  binning  scheme,  as  opposed  to  which  instrument 
data  will  be  used.  Chief  concerns  were  the  following: 

•  The  data  base  should  be  made  as  compact  as  possible.  ITius,  rather  than  including  all  selected 
data  at  instrument  resolution,  the  data  base  should  contain  bin-pass  averages. 

•  The  selection  of  a  spatial  binning  grid  coordinate  system:  Proposals,  presented  at  the  Second 
CRRES-SPACERAD  Science  Team  Meeting  (March,  1990,  San  Diego)  by  Karl  Pfitzer  and  A1 
Vampola.  specify,  respectively,  a  geographic  and  an  L-B  binning  grid,  llie  advantage  of  the 
geographic  grid  is  that  no  magnetic  field  model  is  required  to  define  the  coordinates.  The 
principle  disadvantage  is  that  if  we  later  choose  to  sort  the  geographic  bin-passes  with  respect 
to  an  L-B  or  some  other  coordinate  grid,  an  undesirable  bin  edge  effect  may  appiear  since  a 
geographic  bin  may  straddle  two  or  more  bins  in  the  other  coordinate  system. 

•  The  dimensionality  of  the  snatial  binning:  Should  it  be  one-dimensional,  as  we  used  in  the 
SCATHA  analysis,  or  multi-dimensional?  Although  we  expect  the  radial  variable  (altitude, 
geocentric  distance,  or  L)  to  vary  with  time  much  more  rapidly  than  the  other  coordinates,  a  pass 
through  a  single  radial  bin  may  occasionally  extend  over  a  significant  range  in  one  or  more  of 
the  other  variables  (such  as  latitude  or  local  time).  Thus  a  breakup  of  this  pass  with  respect  to 
the  other  variables  may  be  desired. 

•  Equatorial  or  observed  pitch  angle  binning?  Ihe  data  within  each  spatial  bin  pass  would,  of 
course,  be  further  sorted  and  averaged  within  pitch  angle  bins.  Equatorial  pitch  angle  binning 
has  the  same  objection  as  the  L-binning:  a  magnetic  field  model  is  required  to  transform  from 
the  observed  pitch  angle  to  the  equatorial  pitch  angle. 

1.2. 1.1  Spatial  Binning 

Geographic  bins  are  favored  over  L-B  bins.to  allow  the  user  to  later  select  his  own  variables  with  respect 
to  which  to  sort  the  data.  It  is  much  easier  to  transform  from  geographic  coordinates  to  some  other 
coordinates  than  it  is  to  transform  from  a  B-L  coordinate  system.  To  alleviate  the  bin  edge  effects,  bin- 
widths  would  have  to  be  somewhat  smaller  than  would  be  necessary  if  we  were  binning  directly  with  the 
ultimate  sorting  parameters  (B,  L,  etc.).  The  impact  of  smaller  bin  widths  could  be  mitigated  by  using 
variable  bin  widths. 
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To  aid  in  sizing  the  data  base,  the  number  of  three-dimensional  bins  crossed  by  CRRES  in  one  orbit 
was  determined  for  the  following  scenarios: 


Variable  radial  bin  widths 
Ar  =-  0.01,  1.0  <  r  <  2.0 
Ar  -  0.02,  2.0  <  r  <  6.0 
Ar  *  0.1,  6.0  ^  r  <  10.0 

Fixed  radial  bin  width 

Ar  =  0.01 

The  units  of  r  and  Ar  are  Rg.  In  both  cases  the  geographic  latitude  and  longitude  bin  widths  are  fixed 
at  l®and  15®,  respectively.  The  variable  bin  width  scheme  is  similar  to  A1  Vampola’s  proposal,  except 
that  the  binning  is  in  the  geocentric  radius  r  instead  of  in  L,  and  the  widths  are  5  times  as  large.  Using 
a  model  CRRES  orbit  based  on  preflight  estimates  of  the  initial  conditions,  it  was  found  that  663  bins 
were  crossed  in  the  variable  radial  bin  width  case,  as  opposed  to  1169  in  the  fixed  radial  bin  width  case. 
However,  in  the  variable  bin-width  case,  bins  near  apogee  were  found  to  extend  up  to  40  minutes  in 
duration,  making  identification  and  separation  of  events  from  quiet  periods  difficult.  Furthermore,  local 
time  extents  were  over  one  hour  for  the  large  apogee  bin  widths.  For  these  reasons,  it  seems  prudent 
to  uses  the  smaller  fixed  radial  bin  width. 

As  expected,  the  geocentric  distance  is  the  fastest  changing  variable;  a  3-dimensional  bin  boundaty 
crossing  is  most  often  a  crossing  of  a  radial  bin  boundary.  Although  we  have  chosen  IS®  (~1  hour  in 
local  time)  longitude  bins,  the  longitudinal  interval  within  a  3-dimensional  bin  is  generally  much  less  than 
this,  except  near  perigee.  The  latitudinal  intervals  within  the  3-dimensional  bins  are  considerably  less 
than  1*. 


1.2. 1.2  Sorting  Geographic  Bins  into  L  Bins 

The  geographic  bins  might  later  be  assigned  to  L-bins.  As  mentioned  earlier,  there  could  be  some 
ambiguity  in  such  an  assignment  if  a  geographic  bin  straddles  two  or  more  L-bins.  For  the  b^'o  binning 
cases  mentioned  above,  the  starting  and  ending  L  bin  numbers  of  each  geographic  bin  were  found, 
assuming  that  the  binning  in  L  would  be  analogous  to  the  binning  in  r.  L  was  computed  from  the 
IGRF90  +  Olson-Pfitzer  magnetic  field  model.  It  was  found  that  the  starting  and  ending  L  bins  were 
equal  for  -80  %  of  the  geographic  bins,  thus  permitting  unambiguous  assignment. 


1.2.1.3  Pitch  Ande  Binning 

For  pitch  angle  sorting,  the  observed  pitch  angle  is  preferred  rather  than  the  equatorial  pitch  angle,  since 
the  latter  must  be  computed  from  an  assumed  magnetic  field  model.  For  near-equatorial  latitudes,  the 
two  variables  are  nearly  the  same.  For  higher  latitudes,  a  pitch  angle  bin  of  a  given  width  transforms 
into  a  narrower  equatorial  pitch  angle  bin.  For  example,  in  a  dipole  field,  an  80® -90®  bin  at  20® 
magnetic  latitude  transforms  to  49.3®  -  50.3®  at  the  equator.  Therefore,  the  equatorial  pitch  angle 
resolution  for  particles  mirroring  at  the  observation  point  should  be  excellent. 
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Since  the  angles  a  and  180**  -  a  should  be  equivalent  for  static  modeling,  we  should  sort  these  into  the 
same  bin.  Therefore,  5®  pitch  angle  bins,  centered  at  5®  multiples  from  5®  to  90® ,  are  recommended, 
with  an  additional  loss  cone  bin  0-2.5".  This  would  give  us  19  bins  in  all.  A  single  spin  covers  all  pitch 
angles  from  P  to  180®  -  P  twice,  where  p  is  the  angle  between  the  spin  plane  and  the  magnetic  field 
vector.  If  we  place  supplementary  angles  (a  and  180°  -  a)  in  the  same  bin,  then  all  accessible  bins 
would  be  reached,  at  the  2  rpm  spin  rate,  in  -7.5-15  sec. 


1.2.1.4  Format 

For  each  bin  pass  the  data  base  would  consist  of  the  following: 

Start,  center  and  end  times,  and  the  geocentric  coordinates  at  each; 

For  each  pitch  angle  bin  present; 
the  central  pitch  angle; 
the  total  counts  for  each  selected  channel; 
the  delta  time  for  the  counts  for  each  channel; 
a  data  quality  flag  (contents  TBD). 

In  addition,  for  each  orbit,  there  would  exist  a  small  header  file  indicating  the  start  and  end  times,  the 
selected  instrument  channels,  and  any  other  pertinent  information.  There  would  also  exist  a  separate 
informational  file  containing  generally  constant  information  about  each  of  the  detector  channels  that  may 
be  used  in  the  data  base,  including  its  energy,  and  how  to  obtain  the  flux  from  the  counts.  This  would 
not  be  issued  for  each  orbit,  but  only  a  first  time,  plus  any  other  time  there  are  updates. 

This  format  follows  closely  Karl  Pfitzer’s  proposal.  If  we  use  Al  Vampola’s  suggestion  of  4  energy 
channels  per  decade,  then  we  have  24  channels  in  all,  12  for  protons  and  12  for  electrons.  If  we  assign 
16  bits  to  the  counts  and  16  bits  to  the  delta  time  for  each,  and  we  assign  16  bits  each  to  the  pitch  angle 
and  quality  flag,  then  the  estimated  size,  in  Megabytes  (Mb),  of  the  data  base  for  one  orbit  (variable 
radial  bin  case)  is 

663  X  25  X  19  X  4  X  10-^  1.25  Mb 

The  total  per  week  (16  orbits)  would  be  -20  Mb.  For  the  fixed  radial  bin  specification,  multiply  these 
results  by  1169/663  (=1.76).  This  leads  to  a  weekly  requirement  of  -35  Mb. 


1.2.2  L  vs  Pitch  Angle 

Radiation  belt  populations  have  ainventionally  been  organized  by  their  shell  parameter,  introduced  by 
Mellwain  (1961],  and  commonly  called  L  It  is  often  assumed  that  particles  observed  at  a  given  location 
can  be  assigned  an  L  value  independent  of  pitch  angle.  However,  L  clearly  depends  on  the  mirror  points 
of  a  particle  and,  therefore,  on  its  pitch  angle.  Mcllwain’s  study  of  L  vs  mirror  point  (Figure  15)  shows 
that,  for  an  internal  magnetic  field  model,  the  variation  is  less  than  1%.  However,  CRRES  is 
enanintering  fields  due  to  magnetospheric  currents  not  included  in  the  model  used  by  Mellwain. 
'Hierefore,  a  brief  study  has  been  conducted  of  the  variation  of  L  along  a  field  line  for  the  IGRF85 
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Figure  15.  Deviation  of  L  along  internal  field  model  lines  of  force,  from  Mcllwain  [1961]. 
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(updated  to  17  July  1990)  +  quiet  Olson-Pfitzer  magnetic  field  model. 

L  has  been  calculated  for  field  lines  crossing  the  magnetic  equator  at  rp  =  1.5,  2,  3,  4,  5,  6,  and  7  Rj; 
from  the  offeet  dipole  and  at  local  times  0,  2,  4,  6,  8,  10,  and  12  hr.  The  calculation  was  at  the  mirror 
points  for  equatorial  pitch  angles  5-90°  in  steps  of  5° .  The  orientation  of  the  Earth  was  fixed  such  that 
the  magnetic  longitude  of  the  sun  is  0°.  Figures  16  through  22  show  the  results  for  0  hr  MLT,  and 
Figures  23  through  29  show  the  results  for  12  hr  MLT.  The  parameter  plotted  is 

AL/L^)  =  I(L-Lt)/Lo]  X  100  (%) 

where  Lq  is  the  equatorial  L  value  (the  value  of  L  at  minimum  B).  ITius  is  the  L  value  for  the 
specified  equatorial  pitch  angle.  Each  figure  indicates  at  the  top.  rg.  L^,  MLT,  and  the  minimum  B  value, 
Bq.  For  the  smaller  distances,  the  lowest  pitch  angles  are  excluded  since  the  field  line  trace  stops  at  the 
Earth’s  surface  (the  mirror  points  would,  therefore,  be  below  the  surface). 

Figures  16  and  17  are  consistent  with  Mcllwain’s  results  that  there  is  less  than  1  %  variation.  Recalling 
that  the  Olson-Pfitzer  contribution  is  zero  at  distances  2  Rj^  and  less,  we  are  in  fact  dealing  here  with 
just  the  internal  field.  Figures  18  through  22,  in  which  the  external  field  is  included,  exhibit  a  similar 
shape  which  differ  from  that  of  Figures  16  and  17,  with  maximum  deviations  from  2%  for  r(,  =  3  Rg 
to  just  over  10  %  for  r^  =  7  Rg. 

Figures  23  and  24  (for  12  hr),  like  Figures  16  and  17,  are  at  distances  where  only  the  internal  field 
contributes.  Again  the  variations  are  less  than  1%.  ITie  difference  between  Figures  23  and  24  and  the 
corresponding  midnight  Figures  16  and  17  is  really  due  to  magnetic  longitude  variation,  since  at  these 
distances,  the  magnetic  field  is  independent  of  local  time.  From  3  to  7  Rg  (Figures  25  through  29), 
the  maximum  variation  again  increases  above  1  %  but  is  not  as  high  as  at  midnight.  Unlike  at  midnight, 
the  shapes  of  the  curves  at  these  larger  distances  vary  noticeably  with  distance. 

The  largest  variations  shown  here  are  extreme  cases  which  may  occur  only  rarely  during  the  CRRES 
history  when  the  apogee  is  simultaneously  near  midnight  and  the  equator.  A  preflight  orbital  estimate 
(Figure  30)  indicates  that  for  a  July  19  launch  this  could  occur  about  2  years  after  the  launch,  with  a 
prior  equatorial  crossing  of  apogee  250  days  from  launch  at  20  hours  LT. 


We  have  studied  several  groups  of  models  in  preparation  for  working  with  the  CRRES  data  set.  These 
models  and  the  software  required  to  run  them  constitute  a  whole  suite  of  routines  which  will  be  valuable 
to  CRRES  researchers.  Thus,  we  have  made  them  available  in  a  read  only  library  on  NOS/VE  under 
:CRES.Jordan.CRRES.PLib  for  anyone  to  copy  and  use  at  will.  There  are  five  subdirectories  in  this 
catalog,  grouping  the  software  as  follows: 
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Rgure  16.  Deviation  of  L  along  a  magnetic  field  line  of  force  crossing  the  midnight  equator  at  1.5  Rg.  The  IGRF85 
and  Olson-Pfitzer  quiet  model  field  is  used. 
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Figure  17.  Same  as  Figure  16,  except  for  2  R, 
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Figure  19.  Same  as 
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Figure  20.  Same  as  Figure  16,  except  for  5  Rg. 
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Figure  23.  Same  as  Figure  16,  except  at  noon  instead  of  midnight. 
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Figure  24.  Same  as  Figure  17,  except  at  noon  instead  of  midnight. 
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Figure  26.  Same  as  Figure  19,  except  at  noon  instead  of  midnight. 
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Figure  27.  Same  as  Figure  20,  except  at  noon  instead  of  midnight. 
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Figure  29.  Same  as  Figure  22,  except  at  noon  instead  of  midnight. 
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Figure  30.  Predicted  CRRES  orbital  parameters. 
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1)  Coords  -  coordinate  transformation  routines, 

2)  Eph_Att  --  ephemeris  and  attitude  routines  with  sample  files, 

3)  ExtFlds  -  external  magnetospheric  magnetic  field  models  and  the  code  to  run  them, 

4)  IntFlds  -  internal  magnetospheric  magnetic  field  models  and  the  code  to  run  them,  and 

5)  RadBIts  -  NASA  trapped  radiation  belt  models  and  the  code  to  run  them. 

Within  these  directories,  the  following  routines  reside: 

1)  Coords 

a)  Convrt  --  converts  geodetic  to  geocentric  coordinates  and  vice  versa, 

b)  Transf  --  converts  ECI  to  GSE,  geographic,  GSM.  SM,  VDH,  and  DM,  and  vice  versa, 

2)  Eph_Att 

a)  AGMOD  -  attitude  program  which  performs  line  of  sight  calculations, 

b)  CAF01040  --  sample  attitude  file, 

c)  EOOOSOOl  -  sample  ephemeris  file, 

d)  E000S002  -  sample  ephemeris  file, 

e)  ReadEph  -  reads  an  ephemeris  file  and  prints  out  the  header, 

3)  ExtFlds 

a)  BXYZHV  --  subroutine  containing  the  Hilmer- Voigt  model, 

b)  BXYZMD  -  subroutine  containing  the  Mead  model, 

c)  BXYZNO  --  dummy  subroutine  with  no  external  field  model, 

d)  BXYZOP  --  subroutine  containing  the  Olson-Pfitzer  tilt-dependent  model, 

e)  BXYZT89  -  subroutine  containing  the  Tsyganenko  1989  model, 

0  DrvrEqB  -  sample  driver  program  shows  how  to  set-up  a  model  run;  must  be  linked  to 
a  BXYZ _ file  and  to  the  library  ExMdAft, 

g)  EqBMdJob  -  sample  job  for  submitting  a  model  run  to  the  batch  queue, 

h)  ExMdAft  -  library  of  required  routines  for  a  model  run, 

i)  ExMdSrc  --  source  code  of  the  library  routines, 

j)  HV_In  -  sample  Hilmer- Voigt  input  parameters, 

k)  Md_In  -  sample  Mead  input  parameters, 

l)  T89Par  -  Tsyganenko  1989  parameter  file  required  to  be  in  the  directory  during  a  run, 

m)  T89_In  -  sample  Tsyganenko  1989  input  parameters, 

4)  IntFlds 

a)  CInteg  -  Cain’s  integration  routine  to  obtain  B  and  L, 

b)  GenBL  -  gets  B  and  L  for  a  given  latitude,  east  longitude,  and  altitude  for  one  of  several 

provided  internal  field  models, 

c)  Intrp85  -  Cain’s  interpolation  method  to  get  L  and  the  invariant  latitude  for  1985, 

d)  Intrp90  -  Cain’s  interpolation  method  to  get  L  and  the  invariant  latitude  for  1990, 
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5)  RadBlts 


a)  AESMax 

b)  AESMin 

c)  APSMax 

d)  AP8Min 

e)  Models? 

f)  RadBelt 


-  NASA  trapped  electron  model  for  solar  maximum  conditions. 

-  NASA  trapped  electron  model  for  solar  minimum  conditions, 

-  NASA  trapped  proton  model  for  solar  maximum  conditions, 

-  NASA  trapped  proton  model  for  solar  minimum  conditions, 

“  sample  ''ode  to  run  the  trapped  particle  models,  1987  version,  and 

-  sample  code  to  run  the  trapped  particle  models,  1988  version. 


1.3.2  Other  Developments 

1.3.2.1  Conversion  of  Observed  Pitch  Angles  to  Eouatorial 

The  product-associated  data  base  described  above  provides  distributions  of  trapped  particle  populations 
in  the  observed  pitch  angle  at  the  location  of  the  satellite,  r.  For  further  modeling  and  analysis,  it  is 
useful  to  convert  to  distributions  in  the  equatorial  pitch  angle,  which  would  be  measured  at  the  equatorial 
crossing  of  the  field  line  through  r  (see  Figure  31  below),  thereby  removing  the  variation  of  pitch 
angle  with  latitude.  This  can  be  done  by  the  Liouville  theorem,  from  which  it  follows  that 


/e(«e)  = 


(12) 


where 

/g  =  the  pitch  angle  distribution  at 
=  the  pitch  angle  at  r^; 

/  =  the  pitch  angle  distribution  at  the  point  r, 

a  =  the  pitch  angle  that  a  particle  with  equatorial  pitch  angle  at  would  have  at 

the  point  r. 


From  the  constancy  of  the  first  adiabatic  invariant: 


sin  a  _  B 

sin^Og 


(13) 


where  B  and  B^  are  the  magnetic  field  vector  magnitudes  at  r  and  r^,  respectively. 

Note  that  there  is  a  maximum  pitch  angle  for  which  the  equatorial  distribution  function  can  be 
determined  from  the  distribution  function  /  at  r.  lliis  is  given  by 
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sin^  a. 


C 


(14) 


This  is  the  largest  equatorial  pitch  angle  for  which  particles  can  mirror  at  points  of  latitude  greater  than 
or  equal  to  the  latitude  of  the  point  r  (Figure  31). 

Program  EQMAP  performs  the  conversion,  specified  by  Equation  (13),  of  observed  pitch  angle  to 
equatorial  pitch  angle  for  each  of  19  observed  pitch  angle  bins  for  the  set  of  positional  bin  passes  in  a 
product  associated  data  base.  The  central  and  maximum  angle  of  each  bin  are  converted.  The  bins  are 
assumed  to  be  centered  at  multiples  of  5°,  except  for  a  loss  cone  bin,  0-2..5°,  and  the  last  bin,  87.5° -90°. 
An  input  file  is  required  containing,  for  the  time  center  of  each  positional  bin  pass,  the  magnetic  field 
magnitudes  at  the  satellite  and  at  the  equatorial  crossing.  An  output  file  is  created  containing  all  the 
information  of  the  input  file,  plus  the  resulting  equatorial  pitch  angles. 

1.3.2.1.1  Subroutine  BFLDSM 

Program  OPTRACE  [RADEX,  Inc,  1987)  has  the  capability  to  compute  L  for  a  specified  point.  This 
value  relates  to  the  second  adiabatic  invariant  only  for  particles  which  mirror  at  this  point,  that  is, 
particles  whose  pitch  angles  are  90°  at  this  point.  The  second  adiabatic  invariant  for  a  non-90°  pitch 
angle  particle  is  related  instead  to  the  L  value  at  the  points  where  that  particle  mirrors.  These  points 
are  specified  by 


Bm- 


Bn 


sin^o 


(15) 


where  Bq  is  the  magnetic  field  magnitude  at  our  initial  point,  is  the  mirror  point  magnetic  field 
magnitude,  and  a  is  the  pitch  angle  at  the  initial  point.  ITie  closest  point,  satisfying  Equation  (15)  on 
each  side  of  the  initial  point,  must  be  found.  The  L  value  then  depends  on  the  invariant  integral  between 
these  two  points: 


d 


da 


(16) 


where  c  and  d  are  the  two  points  where  the  particle  mirrors,  B(s)  is  the  magnetic  field  magnitude  at  a 
point  s  on  the  field  line,  is  the  mirror  point  magnetic  field  magnitude  according  to  Equation  (15), 
and  ds  is  an  element  of  length  along  the  field  line. 
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Figure  31.  Magnetic  field  line  passing  through  a  given  vector  location  r  and  crossing  the  Equator  at  r^. 
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To  perform  the  calculation  requires  the  following  steps: 

1)  Trace  the  field  line  through  the  initial  point  to  both  Earth  intercepts; 

2)  Determine  the  B  value,  by  Equation  (15),  for  the  points  where  the  particle  mirrors; 

3)  Search  for  these  points  along  the  field  line  traced  in  step  1; 

4)  Compute  the  invariant  integral  I,  Equation  (16),  for  these  points,  and  then  L  from  Mcllwain’s  or 

Hilton’s  approximation. 

Subroutine  BFLDSM,  a  modification  of  the  OPTRACE  package,  contains  the  option  to  save  a  traced 
field  line  in  a  common  memory  area.  For  each  point  along  the  line,  the  cartesian  coordinates,  altitude, 
and  magnetic  field  magnitude  are  saved.  The  inputs  to  the  routine  are: 

1.  Date  and  GMT 

2.  Positional  coordinates 

The  output  is  the  traced  field  line  just  described,  stored  northern-most  points  first.  Step  2  is  a  simple 
calculation,  using  Equation  (15).  Steps  3  and  4  can  be  performed  by  subroutine  LB,  part  of  the 
BFLDSM  package.  It  takes  as  its  input,  a  traced  field  line,  as  just  described,  and  the  B  value  at  the 
mirror  point  (B^  of  Equation  (15)).  The  output  is  the  desired  L  value.  BFLDSM  has  been  used  in  the 
L  vs  pitch  angle  studies  described  in  section  1.2.2,  above,  and  in  program  PADEPL,  described  below 
(section  1. 3.2.2). 

The  BFLDSM  package  must  be  linked  with  the  OPTRACE  package,  and  with  object  modules  named 
IGRF85  and  BXYZMU,  respectively,  for  the  desired  internal  and  external  field  models.  These  would 
be  used  instead  of  the  default  modules  in  the  OPTRACE  package. 

Incidentally,  the  BFLDSM  package  contains  modules  also  for  the  following: 

1)  Writing  saved  field  lines  to  disk  for  subsequent  display  using  field  line  plot  packages; 

2)  Finding  the  location  and  value  of  the  minimum  magnetic  field  magnitude; 

3)  Finding  the  equatorial  crossing; 

4)  Finding  the  low  altitude  footprints. 
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1.3.2.1.2  CRRES  Coordinate  Transformation  Software  Library 


This  package  permits  one  to  perform  transformations  of  rectangular  coordinates  or  vector  components 
between  the  Earth  Centered  Inertial  (ECI)  coordinate  system  and  the  following: 

Geocentric  Solar  Ecliptic  (GSE) 

Geographic  (GEOG) 

Geocentric  Solar  Magnetospheric  (GSM) 

Solar  Magnetic  (SM) 

Vehicle  Dipole  Horizon  (VDH)  or  Dipole  Meridian  (DM) 

These  coordinate  systems,  described  in  detail  by  Bhavnani  and  Vancour  [1991],  are  summarized  below. 
Since  we  describe  right-handed  coordinate  systems,  it  is  necessary  to  define  only  two  of  the  coordinate 
axes. 


Earth  Centered  Inertial  (ECI) 

The  z-axis  is  along  the  Earth’s  rotational  axis,  positive  north;  the  x-axis  points  to  the  vernal  equinox. 
Origin  is  center  of  the  Earth. 


Geocentric  Solar  Ecliptic  (GSE) 

The  x-axis  points  from  the  center  of  the  Earth  to  the  sun;  the  z-axis  is  perpendicular  to  the  ecliptic  plane 
(parallel  to  the  ecliptic  pole),  positive  north.  Origin  is  center  of  Earth. 

Geographic  (GEOG) 

The  z-axis  is  along  the  Earth’s  rotational  axis,  positive  north;  the  x-axis  points  to  the  Greenwich 
meridian.  Origin  is  center  of  Earth. 

Geocentric  Solar  Magnetospheric  (GSM) 

The  x-axis  points  from  the  Earth  to  the  sun;  the  z-axis  is  perpendicular  to  the  x-axis,  in  the  plane 
containing  the  x-axis  and  the  magnetic  dipole  axis,  positive  north.  Origin  is  center  of  Earth  or  location 
of  offset  dipole. 


Solar  Magnetic  (SM) 

The  z-axis  is  along  the  magnetic  dipole  axis,  positive  north;  the  x-axis  is  perpendicular  to  the  z-axis,  in 
the  plane  containing  the  z-axis  and  the  sun.  positive  toward  the  sun.  Origin  is  center  of  Earth  or 
location  of  offset  dipole. 

Vehicle  Dipole  Horizon  (VDH)  or  Dipole  Meridian  (DM) 

The  z-axis  (H)  is  parallel  to  the  magnetic  dipole  axis,  positive  north;  the  y-axis  (D)  is  perpendicular  to 
the  z-axis  and  the  vector  from  the  origin  to  the  point  of  interest,  positive  east.  Origin  is  point  of  interest. 
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The  Offeet  Dipole  Location 


In  the  GSM,  SM  ,  and  VDH  systems  the  option  exists  to  use  either  the  center  of  the  Earth  or  the 
location  of  the  offeet  dipole  as  the  coordinate  origin.  The  location  of  the  offset  dipole  is  defined  as  that 
position,  which,  if  used  as  a  coordinate  origin,  would  result  in  the  vanishing  of  the  (n=2,m=0)  and 
(n=2,m= 1)  coefficients  of  the  spherical  harmonic  expansion  of  the  Earth’s  magnetic  field  in  a  coordinate 
system  whose  z-axis  is  aligned  with  the  dipole.  (See  Ckapman  and  Bartels  [1940]  for  details).  To  use 
the  offeet  dipole  option,  the  input  vector  (the  vector  in  the  coordinate  system  to  be  transformed  from) 
must  be  in  km.  When  transforming  from  the  ECI  system  to  either  GSM  or  SM,  the  appropriate 
rotation  matrbc  is  applied  to  the  vector  difference  between  the  input  vector  and  the  offeet  dipole 
displacement  from  the  center  of  the  Earth,  in  the  ECI  system.  In  the  inverse  transformation,  the  input 
vector  is  rotated  to  ECI,  then  the  ECI  offeet  dipole  displacement  vector  is  added  to  the  result. 


1. 3.2.2  L  vs  Pitch  Angle  for  Ephemeris 

In  a  static  radiation  belt  model,  the  mirror  point  L  and  B  values  should  serve  as  independent  variables. 
To  develop  such  a  model,  data  should  be  sorted  with  respect  to  these  variables  and  the  energy.  It  is 
frequently  assumed  that  L  is  constant  along  a  field  line,  and  therefore,  independent  of  pitch  angle  for 
particles  observed  at  a  given  position.  However,  as  discussed  previously,  there  can  be  considerable 
variation  near  geosynchronous  altitude  in  a  realistic  magnetic  field  model.  'Fherefore,  it  appears  prudent 
to  obtain  the  pitch-angle-dependent  L  for  each  observation,  rather  than  the  L  at  the  specified  location, 
given  on  the  ephemeris  file.  This  latter  L  is  just  the  value  for  particles  with  90°  pitch  angles.  ITierefore. 
program  PADEPL  has  been  developed  which  provides  a  supplement  to  the  ephemeris  file  containing, 
for  each  record  on  the  ephemeris  file,  the  L  values  for  each  of  the  pitch  angle  bin  centers  defined  in  the 
above  paragraph.  The  pitch-angle  dependent  L  values  may  then  be  computed  for  arbitrary  times  and 
pitch  angles  by  conventional  interpolation  schemes  as  is  done  for  the  standard  ephemeris  variables. 


1.3.2. 3  Adiabatic  Invariants  and  Dipole-Equivalent  Invariants 

To  facilitate  the  study  of  radiation  belt  dynamics.  Schulz  [1991]  (Spring  AGU  Meeting)  has 
recommended  organizing  the  CRRES  radiation  belt  particle  data  in  terms  of  the  adiabatic  invariants. 
Iliese  are  defined  as  action  integrals  pertaining,  respectively,  to  the  basic  three  approximately  separable 
periodic  components  of  the  motion  of  magnetospherically  trapped  particles;  the  spiraling  motion  of  the 
particle  about  a  guiding  center,  the  motion  of  the  guiding  center  along  field  lines,  and  the  longitudinal 
drift  motion  of  the  guiding  center  around  the  Earth.  In  the  absence  of  an  electric  field  component 
parallel  to  the  magnetic  field,  the  first  invariant  reduces  to  a  very  simple  expression: 

p  =  p^/(2mB^).  (17) 

where  p  is  the  particle’s  momentum,  m  its  mass,  and  B^  is  the  mirror  point  magnetic  field  magnitude, 
given  by 


=  B/sin^a, 


(18) 


.W 


where  B  is  the  magnetic  field  magnitude  at  a  given  location  traversed  by  the  particle,  and  a  is  the  pitch 
angle  at  that  location. 

llie  second  adiabatic  invariant  is  given  by 


J  =  2pl„, 


(19) 


where, 


=  /Al-B(s)/B„)ds,  (20) 

with  the  integral  taken  along  the  field  line  between  the  particle’s  mirror  points.  B(s)  is  the  magnetic  field 
magnitude  at  a  point  s  on  the  field  line. 

In  these  expressions,  and  in  what  follows,  the  subscript  m  is  used  to  designate  a  mirror  point  quantity. 
Thus,  is  the  magnetic  field  at  the  mirror  point  of  a  particle,  as  opposed  to  the  point  of  observation. 
I.  without  the  subscript,  can  be  considered  a  function  of  position  if  B^  is  replaced  by  B  at  the  position 
in  question  and  the  integral  is  taken  b)etween  that  position  and  its  conjugate  point.  The  subscript 
emphasizes  the  dependence  of  these  quantities  on  the  mirror  points  rather  than  the  position  of 
observation.  The  mirror  point,  in  turn,  is  a  function  of  the  point  of  observation  and  the  pitch  angle. 

The  third  adiabatic  invariant  is  given  by 


«I»  =  /B  •  da,  (21) 

where  the  integral  is  over  any  surface  enclosed  by  the  drift  orbit. 

Note  that,  if  the  magnetic  field  is  static,  the  particle  momentum  is  fixed  and  the  parameters  and  B^ 
serve  equally  as  well  as  the  adiabatic  invariants  for  parameterizing  the  particle  trajectories,  since  they  are 
then  functions  of  p  and  J  and  are,  therefore,  themselves  invariant.  Thus,  and  B„  specify  the  set  of 
mirror  points  through  which  the  particle  passes  in  its  orbit,  and  therefore,  completely  sp)ecify  the  orbit. 
In  this  case,  the  third  invariant,  the  flux  enclosed  by  the  drift  orbit,  is  not  an  independent  parameter. 
However,  its  place  as  a  constant  of  the  motion  may  be  taken  by  the  particle’s  energy.  The  constancy  of 
these  parameters  disappears  in  time-varying  fields  in  which  the  resulting  induction  electric  field  does  work 
on  the  particles,  thus,  changing  their  energies  and  momenta.  and  B^  must  correspondingly  be 
adjusted  to  conserve  p  and  J. 

McRwain  [1961]  has  defined  a  shell  parameter  as  a  function  of  and  B^.  This  parameter  is, 
therefore,  adiabatically  conserved  if  and  B„  are.  In  a  dipole  field,  is  the  equatorial  crossing 
distance  of  the  shell  of  field  lines  through  which  the  particle  passes.  Therefore,  at  a  given  location,  all 
particles  have  the  same  value  regardless  of  pitch  angle.  Since  this  relationship  remains  approximately 
valid  in  the  Earth’s  internal  field,  has  a  distinct  advantage  over  I^  for  organizing  the  inner  belts. 
As  a  function  of  position,  L  is  constant  along  a  dipole  field  line  and  approximately  constant  along  field 
lines  of  the  Earth’s  internal  field. 
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1.3.2.3.1 


A  problem  with  the  adiabatic  invariants  is  the  difficulty  in  working  with  them  due  to  lack  of  familiarity. 
We  are  more  accustomed  to  energy  spectra,  radial  distributions,  and  pitch  angle  distributions,  llius,  we 
are  motivated  to  define  a  set  of  equivalent  invariants  -  the  energy,  L„^nd  equatorial  pitch  angle  that 
particles  with  a  given  set  of  adiabatic  invariants  would  attain  if  the  magnetic  field  changed  adiabatically 
to  a  specified  reference  field.  In  the  past,  the  Mead  model  has  been  used  as  the  reference  field  [Boas, 
et.  al.,  1989].  A  simpler  choice  would  be  a  dipole  field.  We  would,  therefore,  call  these  parameters 
dipole-equivalent  invariants.  First,  we  have,  for  the  dipole-equivalent 

Lmd  =  (22) 

where  M  is  the  dipole  moment  of  the  chosen  dipole  reference  field.  A  review  of  Roederer  [1970]  and 
Schulz  and  Lanzerotti  [1974]  leads  to  the  conclusion  that  this  is  the  L  parameter  used  by  them  in 
diffusion  equations,  rather  than  Mcllwain’s  parameter.  Unlike  the  latter,  is  a  function  of  the  third 
adiabatic  invariant,  and,  therefore,  would  not  change  under  adiabatic  variations  of  the  magnetic  field. 
Furthermore,  this  is  actually  independent  of  a  particle’s  mirror  points  in  a  symmetric  non-dipole 
field.  Thus,  any  dependence  on  mirror  points  reflects  shell-splitting  in  an  asymmetric  field. 

To  obtain  the  dipole-equivalent  equatorial  pitch  angle,  we  note  the  following  [Roederer,  1970]: 

=  constant  (23) 

which  follows  from  Equations  (17)  and  (18),  and  constancy  of  J  and  p.  After  substituting  known  dipole 
expressions  for  and  as  functions  and  a^,  we  solve  for  Oj.  Finally,  constancy  of  the  first 
adiabatic  invariant  supplies  the  ratio  p^/p,  where  p^  is  the  dipole-equivalent  momentum,  and  p  is  the 
observed  momentum.  From  this  ratio,  we  may  obtain  the  dipole-equivalent  energy  for  any  observed 
energy. 

1.3.2.3.2  Software 

Figure  32  summarizes  software  that  has  been  developed  to  compute  the  adiabatic  invariants  and  dipole- 
equivalent  invariants.  As  indicated  in  the  top  block,  we  are  given  the  particle  position,  pitch  angle, 
energy,  mass,  and  the  time.  The  first  step  is  to  compute  the  set  of  energy-independent  parameters  for 
specified  positions  and  pitch  angles  in  a  given  field  model.  These  are  the  quantities  1^^^,  L^,  and  ♦. 
Subroutine  ADINV  performs  this  function  (in  this  and  all  subsequent  blocks,  CAPS  enclosed  by 
parentheses  indicate  subroutine  names).  This  is  the  most  CP-time-intensive  step  of  the  process,  and  will 
be  described  in  more  detail  below.  The  remainder  of  the  procedure  has  two  alternatives.  If  one  desires 
the  adiabatic  invariants,  it  is  prudent  to  take  the  left-hand  branch,  in  which  the  first  two  adiabatic 
invariants,  J  and  p.  are  obtained  by  direct  substitution  (subroutine  ADINV2)  of  and  B^  into 
Equations  (17)  and  (19).  The  dipole-equivalent  invariants  are  then  obtained,  as  described  above,  using 
subroutine  DIPINV.  The  right-hand  branch  may  be  followed  if  one  needs  only  the  dipole-equivalent 
invariants,  but  is  computing  them  for  several  energies  and/or  masses.  Then  it  is  noted  that  the  first  two 
of  these,  and  Oj,  are  independent  of  mass  and  energy,  thus,  needing  to  be  determined  only  once  at 
a  given  location,  time,  and  pitch  angle.  Along  with  these,  a  third  mass/energy  independent  parameter, 
the  ratio  p/pj  is  also  determined.  ITiese  are  collectively  called  energy-independent  diptrle-equivalent 
parameters.  Finally,  the  dipole-equivalent  energy  can  be  determined  for  each  mass  and  energy 
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Figure  32.  Overview  of  software  for  cpmputation  of  adiabatic  invariants  and  dipole  equivalent  invariants. 
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(subroutine  DIPINVB). 

The  computation  of  the  energy-independent  parameters  (subroutine  ADINV)  requires  three  stages: 

1)  Trace  the  magnetic  field  line  through  the  position  of  interest  and  compute  the  magnetic  field  value 

at  that  specified  location. 

2)  For  each  pitch  angle,  determine  the  mirror  point  magnetic  field  value  and  the  mirror  point 

integral  using  Equations  (18)  and  (20). 

3)  For  each  pitch  angle,  compute  the  third  invariant  ♦  Equation  (21). 

Stage  3  requires  the  determination  of  the  drift  shell,  defined  by  the  set  of  field  lines  containing  the  mirror 
points  of  the  requisite  and  B^^.  llie  magnetic  equatorial  crossing  distances  of  the  drift  shell  are 
determined  at  each  of  25  equally  spaced  local  times  beU\’een  0  and  12  hours.  That  portion  of  the  shell 
between  12  and  24  hours  is  determined  from  the  portion  between  0  and  12  hours  by  dawn-dusk 
symmetry.  The  required  field  line  crossing  the  equator,  at  a  given  local  time,  is  found  by  iterative  search 
(subroutine  LBSRCH).  Ihe  initial  guess  for  the  crossing  distance  is  the  value  of  Mcllwain’s  function 
of  B|^  and  I^.  Then  the  value  of  L  at  the  point  B=B|jj  along  this  field  line  is  found.  If  this  value  is 
within  10'^  of  the  required  L^,  then  the  search  is  terminated.  Otherwise,  the  difference  is  used  to 
correct  the  equatorial  crossing  distance  and  the  process  repeated  until  the  crossing  distance  is  determined 
to  within  specified  accuracy  (10"^  R|;).  Generally  this  takes  3-5  iterations. 

Given  the  set  of  equatorial  crossing  distances  defining  the  drift  shell,  the  third  invariant  is  computed  by 
(subroutine  INVAR3),  using  code  provided  by  Lockheed,  Inc.  {Nightingale,  et.  al.,  1986).  ITiis  code 
computes  the  third  invariant  in  3  different  ways: 

1)  Compute  directly  the  magnetic  flux  crossing  that  portion  of  the  equatorial  surface  bounded  by  the 

drift  shell.  In  this  approach,  the  region  of  the  equatorial  surface  near  the  origin  is  replaced  by  a 
hemisphere  to  avoid  the  singularity  in  the  magnetic  field  at  the  origin. 

2)  Compute  the  flux  crossing  that  portion  <  f  a  hemispherical  surface  of  specified  radius  enclosed  by  the 

drift  shell. 

3)  Compute  the  line  integral  of  the  vector  potential  along  the  intersection  of  the  drift  shell  with  a 

hemispherical  surface  of  specified  radius. 

For  a  divergence-free  magnetic  field,  these  approaches  give  the  same  result.  However,  since  the  latter 
two  methods  require  the  additional  field  line  tracing  to  the  specified  surface,  only  the  first  method  has 
been  used  here. 


54 


1.4  References 


Baker,  D.  N.,  Higbie,  P.  R..  Hones,  E.  W.,  Jr.,  and  Belian,  R.  D.,  High  resolution  energetic  particle 
measurements  at  6.6  R^:  3.  Low-energy  electron  anisotropies  and  short-term  substorm  predictions,  J. 
Geophys.  Res.,  vol.  83,  no.  AlO,  p.  4863,  1978. 

Bass,  J.  N.,  Gussenhoven,  M.  S.,  and  Redus,  R.  H..  "The  Importance  of  Adiabatic  Variations  in  Trapped 
Particle  Distributions  Observed  by  the  SCATHA  Satellite",  GL-TR -89-0233,  29  Nov  1989.  ADA220^9. 

Bhavnani,  K.  H.,  and  Vancour,  R.  P.,  "Qwrdinate  Systems  for  Space  and  Geophysical  Applications",  PL- 
TR-91-22%,  Dec  1991. 

Chapman.  S.,  and  Bartels,  J.,  Geomagnetism,  Vol.  II,  Oxford  at  the  Qarendon  Press.  London,  p.  648, 
1940. 

Gussenhoven,  M.  S.,  Mullen,  E.  G.,  and  Sagalyn,  R.  C,  "CRRES/SPACERAD  Experiment 
Descriptions",  Air  Force  Geophysics  Laboratory  Report  AFGL-TR -85-0017,  p.  66,  1985,  ADA160504. 

Hanser,  F.  A,  Hardy.  D.  A,  and  Sellers,  B..  “Calibration  of  the  Rapid  Scan  Particle  Detector  Mounted 
in  the  SCATHA  Satellite",  Air  Force  Geophysics  Laboratory  Report  AFGL-TR-79-0167,  1979, 
ADA082382. 

Heinemann,  M.,  Representations  of  currents  and  magnetic  fields  in  anisotropic  magnetohydrostatic 
plasma,  J.  Geophys.  Res.,  vol.  95,  no.  A6.  p.  7789, 1990. 

Hilmer,  Robert  V.,  A  Magnetospherie  Magnetic  Field  Model  with  Flexible  Internal  Current  Systems, 
A  thesis  submitted  in  partial  fulfillment  of  the  requirements  for  the  degree  doctor  of  philosophy.  Rice 
University,  Dept,  of  Space  Physics  and  Astronomy,  Houston,  TX,  1989. 

Mcllwain,  C  E.,  "Coordinates  for  Mapping  the  Distribution  of  Magnetically  Trapped  Particles",  J. 
Geophys.  Res.,  Vol.  66,  No.  11,  Pp.  3681-3691.  1961. 

Mullen,  E,  G..  and  Gussenhoven,  M.  S.,  "SCAITIA  Environmental  Atlas",  Air  Force  Geophysics 
Laboratory  Report  AFGL-TR -83-0002,  pp.  27-76,  1983,  ADA131456. 

Nightingale,  R.  W.,  Chiu,  Y.  T.,  Davidson,  G.  1'..  Francis,  W.  E.,  Rinaldi,  M.  A,  Robinson,  R.  M.,  and 
Vondrak.  R.  R.,"A  Space  Radiation  Test  Model  Study",  AFGL-'l’R -86-0064,  1986,  ADA169412. 

Olson,  W.  P.,  and  Pfitzer.  K.  A,  "Magnetospherie  Magnetic  Field  Modeling",  Annual  Scientific  Report, 
AFOSR  Contract  No.  F44620-75-C-0033,  McDonnell  Douglas  Astronautics  Company,  Huntington 
Beach.  CA  1977. 

Radex,  Inc,  "CRRES  Data  Processing  Ephemeris  File  Generation  System  Product  Development 
Specification".  Radex,  Inc.,  Bedford,  MA  01730,  RXR87041,  1987. 


55 


Radex.  Inc,  "Integrated  Analytical  Evaluation  and  Optimization  of  Model  Parameters  Against 
Preprocessed  Measurement  Data”,  Radex,  Inc  Final  Report,  GL-TR-89-0163,  1989,  ADA215518. 

Roederer,  J.  G.,  "Dynamics  of  Geomagnetically  Trapped  Radiation",  Springer-Verlag,  Berlin/Heidelber ’y 
New  York.  1970. 

Schulz,  M.,  and  Lanzerotti  .L  J..  "Particle  Diffusion  in  the  Radiation  Belts",  Springer-Verlag, 
Berlin/Heidelberg/New  York,  1974. 

Schulz.  M.,  and  McNab.  M.  C,  "Source-Surface  Model  of  the  Magnetosphere",  Aerospace  Corp.,  Space 
Sciences  Lab,  Rpt.  No.  SSL-87(7165)-1,  Nov  1986. 

Schulz.  M..  "Opportunities  for  Radiation  Belt  Modeling:  Empirical  and  Theoretical”,  Spring  Meeting 
1991  Program  and  Abstracts,  EOS  Suppl.,  p.  232,  1991. 

Stern,  D.  P.,  "Parabolic  Harmonics  in  Magnetospheric  Modeling:  ITie  Main  Dipole  and  the  Ring 
Current".  J.  Geophys.  Res.,  Vol.  90.  No.  All,  Pp.  10.851-10,863.  1985. 

Tsyganenko,  N.  A..  A  magnetospheric  magnetic  field  model  with  a  warped  tail  current  sheet.  Planet. 
Space  Sci.,\'o\.  37,  no.  1,  p.5,  1989. 

Vasyliunas,  V.  M.,  Mathematical  models  of  magnetospheric  convection  and  its  coupling  to  the 
ionosphere,  in  Particles  and  Fields  in  the  Magnetosphere,  edited  by  B.  M.  McCormac.  p.  60,  D.  Reidel, 
Hingham.  Mass.,  1970. 

We.st.  H.  L..  Jr..  The  signatures  of  the  various  regions  of  the  outer  magnetosphere  in  the  pitch  angle 
distributions  of  energetic  particles,  in  Quantitative  Modeling  of  Magnetospheric  Processes,  Geophys, 
Monogr.  Ser..  vol.  21,  edited  by  W.  P.  Olson,  p.  150.  AGU,  Washington,  D.  C,  1979. 


56 


2.0  DYNAMIC  MAGNETIC  FIKLD  MODKLii 


.?.l  Models  Published  Before  198^ 


In  ihe  Earth’s  magnetosphere,  the  magnetic  field  is  the  sum  of  two  parts:  that  which  we  will  aill  the 
"internal”  field  because  it  is  due  to  sources  (currents)  in  the  Earth’s  interior,  and,  that  which  we  will  call 
the  "external"  magnetic  field  because  it  is  due  to  sources  outside  the  Earth's  surface.  Since  its  sources  are 
in  the  interior  of  the  Earth,  the  internal  field,  at  a  point  outside  the  Earth’s  surface,  is  curl-free. 
I'herefore,  the  internal  field  can  be  expressed  as  the  gradient  of  a  scalar  potential. 

Since  the  magnetic  field  is  also  divergence-free,  the  scalar  potential  satisfies  Laplace’s  equation  and  is, 
therefore,  commonly  expressed  as  an  expansion  in  spherical  harmonics. 

The  external  field,  unlike  the  internal,  is  not,  in  general,  current-free  and  usually  cannot  be  expressed  as 
the  gradient  of  a  scalar  potential.  However,  since  the  external  currents  are  more  amenable  to 
observation,  more  can  be  said  about  them  than  the  internal  currents.  Ilius,  external  field  modeling  has 
frequently  begun  with  models  of  the  current  from  which  the  field  is  subsequently  derived.  The  currents, 
included  in  most  models,  are  those  known  as  ring  currents,  cross-tail  currents,  and  magnetopause 
currents.  ITie  ring  currents,  carried  by  trapped  ions  and  electrons  in  the  20-.S0  keV  energy  range, 
circulate  longitudinally  around  the  Earth  at  distances  2-7  R^;.  The  cross-tail  currents  flow'  from  dawn 
to  dusk  across  the  night-side  plasma  sheet.  Some  models  also  include  the  associated  return  currents, 
flowing  around  the  magnetopause  from  dusk  to  dawn.  The  magnetopause  currents  are  the 
Giapman-Ferraro  boundary  currents,  which  are  formed  by  the  interaction  of  the  magnetosphere  with 
the  solar  wind.  Ihe  field,  due  to  the  magnetopause  currents  alone,  can  be  expressed  in  terms  of  a  scalar 
potential  (for  points  inside  the  magnetosphere).  Ihere  is  another  system  of  currents,  the  field-aligned 
currents,  which  none  of  the  models  to  be  discussed  here  include  explicitly. 

Ihe  CRRES  satellite  will  provide  data  for  a  generation  of  new  models  of  the  Earth's  radiation  belts. 
Since  the  orbits  of  the  trapped  high  energy  radiation  belt  particles  are  controlled  by  the  Earth's  magnetic 
field,  the  accuracy  of  the  models  will  largely  l>e  determined  by  the  accuracy  of  the  magnetic  field  models 
used  in  ordering  the  data.  Present  radiation  belt  models  neglect  the  external  magnetic  field.  While  this 
approximation  may  be  acceptable  for  the  inner  belt,  it  becomes  questionable  at  larger  distances  beatuse 
the  external  field  accounts  for  an  increasing  fraction  of  the  total  field. 

Ihis  subsection  reviews  the  models  listed  in  Table  2.  We  begin  with  the  derivation  and  mathematical 
construction  of  these  models.  A  aimparison  of  computation  times  for  the  models  is  given. 

Table  2.  External  Magnetic  Field  Models 

Olson-Pfitzer  Tilt  Dependent  (Olson  and  Pfitzer.  1V77| 

Olson-Pfitzer  Dynamic  (Pfitzer,  et.  al..  1988} 

Mead-Fairfield  (Mead  and  Fairfield.  197.5 j 
Tsyganenko-Usmanov  [  Taj/janento  and  Usmanov.  1982| 

Tsyganenko  (Tsyganenko.  1987| 

Stern  Parabolic  Magnetopause/Strelched  Magnetosphere  |  Stern,  198.5,  198'7) 


57 


Primary  interest  for  radiation  belt  analysis  are  the  properties  of  the  Mcllwain  L  parameter  and  the 
magnitude  of  the  total  (internal  +  external)  magnetic  field  vector  since  these  have  been  the  traditional 
sorting  parameters  for  modeling  the  trappied  particle  populations.  Comparisons  of  the  L  parameter 
among  the  models  are  presented.  Next,  we  compare  the  total  magnetic  field  vector  derived  from  the 
models  with  those  observed  by  the  SCATHA  satellite  for  a  quiet  and  a  moderately  active  period.  Finally, 
model  depiction  of  the  night  side  equatorial  magnetic  field  depression  as  a  function  of  distance  is 
compared  with  averaged  measurements  by  the  AMP'fE  and  OGO  satellites.  The  equatorial  field  is 
critical  to  the  magnetic  mapping  of  auroral  phenomena  to  the  plasma  sheet. 


2.1.1  Model  Derivation  and  Description 


Most  of  the  models  begin  with  some  description  of  the  three  current  systems  we  mentioned  above,  viz., 
the  ring  currents,  the  cross-tail  currents,  and  the  magnetopause  currents.  The  contributions  of  these 
three  current  systems  are  modeled  separately  and  summed  to  produce  the  total  magnetic  field  vector. 

In  this  report,  we  will  frequently  refer  to  the  Solar  Magnetic  (SM)  and  Geocentric  Solar  Magnetospheric 
(GSM)  coordinate  systems,  following  the  terminology  of  Russell  [1971].  In  both  systems  the  x-z  plane 
contains  the  dipole  and  the  Earth-Sun  line,  while  the  y-axis  is  positive  toward  dusk.  In  the  SM  system 
the  z-axis  is  parallel  to  the  dipole,  positive  toward  the  north,  while  the  x-axis  is  positive  on  the  sunward 
side  of  the  dipole.  In  the  GSM  system,  the  x-axis  is  along  the  Earth-Sun  line,  positive  toward  the  Sun. 
while  the  z-axis  is  positive  on  the  northward  side  of  the  Earth-Sun  line. 


2. 1.1.1  Olson-Pfitzer  Tilt  Dependent  Model 

Ilie  Olson-Pfitzer  tilt  dependent  model  (O/son  and  Pfitzer.  1977]  is  a  revision  of  the  tilt-averaged  model 
of  Olson  and  Pfitzer  [l974j.  l  ilt  here  refers  to  the  aimplement  of  the  Sun-dipole  angle.  ITe 
motivation  for  the  revision  was  that  the  tilt-averaged  model  is  not  a  good  representation  of  the  zero  tilt 
case,  primarily  at  large  distances  on  the  midnight  equator.  When  the  dipole  is  tilted,  the  plasma  sheet 
has  l^een  observed  to  be  parallel  with  the  GSM  equatorial  (x-y)  plane  but  not  coincident  with  it  (Stern, 
1976],  Simple  models  \\furayama.  1966;  Speiser  and  Ness.  1%7;  Fairfield  and  Ness.  1970]  postulate 
that  the  plasma  sheet  is  hinged  to  the  dipole  equatorial  plane  8-11  from  the  Earth.  Locations  on 
the  dipole  equator  outside  the  hinging  distance  are.  therefore,  in  the  plasma  sheet  only  when  the  dipole 
is  imtilted.  Since  the  plasma  sheet  is  an  area  of  minimum  field  strength,  the  field  at  zero  tilt  is  generally 
smaller  than  the  average  field  for  such  points  on  the  dipole  equator. 

Ifie  SM  components  of  are  expressed  as  sums  of  terms  of  the  form 

*  fljjkexp(  0.0tir^)|z'y'z'‘ 

where  the  arefficients  and  8-^  are  simple  linear  or  quadratic  functions  of  the  tilt,  llie  model 
contains  no  dependence  on  magnetic  activity  but  is  proposed  as  an  average  model  for  quiet  conditions 
(Kp  =  0.  1).  Spatially,  it  is  valid  from  2.5  Rj;  to  1.5  R^;,  except  beyond  the  noon  magnetopause  at  10.5 
R(_.  llte  model  is  not  rigorously  divergence  free. 


2,1. 1.1.1  Model  Derivation 


Like  the  1974  tilt-averaged  model,  the  tilt-dependent  model  was  constructed  by  first  considering  models 
for  the  current  systems.  'Fhe  magnetopause  currents  are  similar  to  those  developed  by  OUon  [1969]  from 
a  self-consistent  solution  to  the  Mead-Beard  pressure  balance  condition.  However,  for  the  tilt-dependent 
model,  an  empirical  shape  for  the  magnetopause  surface  was  chosen  instead  of  Olson’s  self-consistent 
solution.  Magnetopause  currents  and  their  resulting  contribution  to  the  field  were  computed  numerically. 

ITie  ring  and  cross-tail  currents  are  constructed  from  a  system  of  wires  [Olson,  1974].  The  ring  currents 
were  made  from  appropriately  placed  elliptical  wires,  such  that  their  noon  crossings  were  at  lower  radial 
distances  than  their  midnight  crossings,  llierefore,  the  resulting  ring  current  is  not  azimuthally 
symmetric.  However,  dawn-dusk  symmetry  was  retained  by  placing  the  elliptical  foci  on  the  SM  x-axis. 
The  flow  in  the  inner  three  wires  is  eastward,  while  in  the  remaining  wires,  it  is  westward.  The  boundary 
between  eastward  and  westward  current  density  is  ~5  Rg  at  midnight,  3  at  noon.  The  cross-tail 
currents  are  constructed  from  loops  consisting  of  dawn-to-dusk  flow  on  the  northern  and  southern 
surfaces  of  the  hinged  plasma  sheet  =  ±3  Rg),  coupled,  respectively,  to  return  currents  over  the 

northern  and  southern  surfaces  of  the  magnetopause.  The  loops  are  tilted  relative  to  the  normal  to  the 
plasma  sheet  in  order  to  model  the  decay  of  with  increasing  distance  down  the  tail.  When  the  dipole 
is  tilted,  the  loops  near  the  inner  edge  are  tilted  toward  the  normal  to  the  dipole  equator  plane. 

Parameters  specifying  the  wires  (position,  orientation,  size,  and  shape),  their  currents,  and  the 
magnetopause  shape  were  adjusted  so  that  the  resulting  total  external  field  fit  OGO-3  and  -5  observations 
in  the  inner  magnetosphere  [5«yiura,  et.  al.,  1971],  and  Explorer-33  and  -35  observations  in  the  tail 
[Mihalov,  et.  al.,  1968;  Behannon,  1970;  Meng  and  Mihalov,  1972].  The  resulting  field,  still  numerical 
in  form,  was  fit  by  linear  least  squares  to  the  analytic  series  described  above. 


2.1. 1.2  Olson-Pfitzer  Dynamic  Model 

This  model  [Pfitzer,  et.  al.,  1988]  is  the  result  of  a  series  of  event  studies  performed  for  NASA 
Coordinated  Data  Analysis  Workshops  (CDAW).  'llie  principal  feature  of  this  model  is  the  introduction 
of  variable  strength  factors  multiplying  the  quiet  models  of  the  fields  of  the  three  current  systems.  The 
magnetopause  strength  factor  varies  inversely  as  the  cube  of  the  standoff  distance,  which  may  be 
computed  from  the  solar  wind  pressure  and  speed  by: 

Rj,  =  98/(nv2)>^  (24) 

with  Rjj  in  Rg,  n  the  solar  wind  ion  density  in  particles/cc.  and  v  the  solar  wind  speed  in  km/sec.  The 
quiet  model  assumes  =  10.5  Rg.  The  field  due  to  the  magnetopause  is  scaled  geometrically  in 
accordance  to  this  scaling  of  the  magnetopause.  'Hius.  the  field  at  the  position  r  due  to  the  scaled 
magnetopause  is 
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where  k  =  R^^/10.5. 


(25) 


The  contribution  due  to  the  tail  current  system  is  scaled  in  the  same  way  as  the  magnetopause 
contribution,  both  in  strength  and  geometry.  An  index  suggested  by  Akosofii  I  Olson  and  Pfitzer,  1932], 
which  was  successfully  employed  in  tail  field  modeling  for  the  CDAW  2  event,  was  found  to  be 
unsatisfactory  for  subsequent  events.  The  ring  current  system  is  scaled  in  strength  only,  using  the 
index: 


^ring 


=  1.0  -  0.03  a. 


(26) 


Ilie  unsealed  (quiet)  models  of  these  systems  are  similar  in  mathematical  form  to  the  quiet  tilt- 
dependent  model.  However,  the  dynamic  mode!  is  valid  only  for  zero  tilt.  Spatially  the  model  is  valid 
from  2.5  Rg  to  60  Rg  on  the  night  side,  and  out  to  the  magnetopause  on  the  day  side.  Like  the  tilt- 
dependent  model,  the  dynamic  model  is  not  divergence-free. 


2.1.1.3  Mead-Fairfield  Model 

In  the  Mead- Fairfield  model  [Mead  and  Fairfield,  1975]  each  component  of  B  is  expressed  as  an 
expansion  in  products  of  polynomials  in  the  SM  coordinates  and  the  dipole  tilt,  llie  SM  coordinates 
are  actually  rotated  westward  by  the  approximate  solar  wind  aberration  angle  of  4°  (the  shift  in  the  solar 
wind  direction  due  to  the  orbiting  Earth).  'Hius.  there  is  symmetry  about  a  plane  4°  west  of  the 
noon-midnight  meridian,  rather  than  about  the  noon-midnight  meridian  itself.  For  each  of  4  bins 
(Q.C*";  0-2';  2-9;  3-9)  the  coefficients  were  derived  from  the  data  by  least  squares,  enforcing  the  zero 
divergence  condition  by  a  set  of  linear  constraints  on  the  coefficients.  The  model  provides  no  description 
of  the  individual  current  systems  but  the  total  current  can  be  derived  by  taking  the  curl  of  the  magnetic 
field  vector. 

'Hie  data  base  consisted  of  12,616  averaged  vector  measurements  by  magnetometers  on  board  satellites 
Explorer  33  (IMP  A-D).  Explorer  34  (IMP-4),  Explorer  41  (IMP-5),  and  Explorer  43  (lMP-6) 
[Behannon,  1968;  Fairfield,  1%9;  Fairfield  and  Ness,  1972;  Fairfield,  1974],  Hie  averaging  was  over 
half-Earth  radii,  which  typically  spanned  10-15  minutes.  'Hie  data  base  contains  gaps  in  the  near-Earth 
equatorial  region,  and  at  high  southern  latitudes.  Specifically,  there  are  few  measurements  within  5  Rg 
of  the  Earth  and  none  within  4  Rg.  Ihus,  the  model  is  valid  for  radial  distances  between  5  and  15  Rg. 


2.1. 1.4  Tsyganenko-Usmanov  and  Tsvganenko  Models 

The  Tsyganenko-Usmanov  [Tsyganenko  and  Usmanov,  1982]  and  Tsyganenko  [1987]  models  are 
expressed  as  sums  of  the  contributions  from  the  three  external  sources  mentioned  above.  Tlie  ring  current 
and  tail  current  contributions  are  expressed  analytically  by  functions  which  are  the  curls  of  vector 
potentials;  hence,  their  divergences  vanish  identically.  Both  models  employ  a  mathematically  simple 
azimuthally  symmetric  ring  current  circulating  around  the  dipole.  In  contrast  to  the  Olson-Pfitzer 
models,  the  ring  current  flow  direction  is  westward,  independent  of  distance  from  the  Earth.  Accurate 
modeling  of  the  innermost  region  was  not  attempted  due  to  insufficient  data  coverage.  In  the 
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Tsyganenko-Usmanov  model,  the  tail  current  flows  mainly  from  dawn  to  dusk  on  a  flat  sheet  "hinged" 
or  offeet  a  distance  on  the  GSM  z-axis  due  to  the  dipole  tilt.  Both  models  use  tilt-dependent  products 
of  polynomials  and  decaying  exponentials  in  the  coordinates  to  represent  the  magnetopause  current 
contribution.  For  these,  the  zero  divergence  condition  was  enforced  while  determining  the  coefficients. 
The  data  bases  for  these  models  contain  the  Mead-Fairfield  data  base  as  a  subset.  The  Tsyganenko  data 
base  is  more  extensive  at  large  distances  down  the  tail;  hence,  it  has  a  more  elaborate  tail  model. 

2.1. 1.4.1  Ring  Current  Model 

The  ring  current  field  in  both  models  is  given  by: 


flp  =  flg[12pV/(p''  +  ^  4)5^] 

(27a) 

=  4flg((2z''  -  p'"  +  8)/(p''  +  ^  4)5^] 

(27b) 

o 

II 

0? 

(27c) 

Where  p'  =  p/pg,  z'  =  z/pg,  p,  z,  and  <J)  are  cylindrical  coordinates  with  the  z-axis  parallel  to  the 
dipole,  and  Bg  and  pg  are  adjustable  parameters.  The  maximum  ring  current  density  is  located  at 
p  =  Pg,  while  the  strength  of  the  current  is  proportional  to  Bg. 

2.1. 1.4.2  Tail  Current  Models 

The  tail  current  field  is  the  product  of  a  field  due  to  an  infinitely  extended  current  in  the  dawn-dusk 
direction  and  a  y-dependent  factor  that  was  introduced  to  improve  model-observation  agreement  in  the 
dawn  and  dusk  regions  of  the  near-Earth  magnetotail.  The  Tsyganenko-Usmanov  model  contains  seven 
adjustable  parameters  and  the  Tsyganenko  model  contains  nine.  The  Tsyganenko  model  also  contains 
return  current  sheets  at  30  above  and  below  the  GSM  equator,  each  equal  in  magnitude  to  half  the 
main  tail  currents. 

2.1. 1.4.3  Magnetopause  Current  Models 

The  magnetopause  current  expansion  consists  of  products  of  powers  of  the  coordinates  y  and  z 
perpendicular  to  the  Earth-Sun  line,  decaying  exponentials  down  tail,  and  the  sine  and  cosine  of  the  tilt 
angle.  There  is  no  explicit  physical  connection  of  the  parameters  with  any  current  structure;  thus,  the 
zero  divergence  criterion  is  not  automatically  satisfied  for  any  set  of  chosen  parameter  values.  Thus,  the 
divergence  requirement  was  imposed  by  a  set  of  linear  constraints  on  the  parameters  during  the  least 
squares  fitting  procedure.  There  are  18  linear  coefficients  in  the  Tsyganenko-Usmanov  magnetopause 
model  and  a  single  nonlinear  parameter  specifying  the  exponential  decay  down  the  tail.  The  Tsyganenko 
model  contains  22  linear  parameters  and  two  nonlinear  parameters.  The  authors  assert  that  these  terms 
also  include  average  contributions  of  field  aligned  currents. 
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2. 1.1. 4.4 


The  Tsyganenko-Usmanov  data  base  consists  of  the  Mead-Fairfield  data  set  and  data  from  the  HEOS-1 
and  HEOS-2  spacecraft  \Htdgtcock  and  Thomas.  1975].  ITie  additional  data  are  6248  vector  averages 
in  the  radial  range  6-35  Rg.  In  the  Tsyganenko  data  base,  this  is  augmented  by  11,150  similar  averages 
from  six  IMP  spacecraft  (A,  C,  E,  F,  G,  and  I)  in  the  range  -66  Rg  <  <  -15  Rg,  and  6675  averages 

from  IMP-H  and  IMP-J  at  down-tail  distances  25-45  Rg.  Apparently,  this  additional  down  tail  data 
motivated  the  refined  Tsyganenko  tail  model.  Of  particular  note  is  that  the  additional  data  (in  both 
models)  do  not  fill  in  the  near-Earth  gap  left  by  the  Mead-Fairfield  model  data  base.  Therefore,  these 
two  models  must  also  be  suspect  within  5  Rg  of  the  Earth.  The  minimal  coverage  of  the 
Tsyganenko-Usmanov  data  set  beyond  20  Rg  down  tail  defines  the  limit  for  that  model,  llie 
Tsyganenko  model  is  supplied  in  two  versions:  a  long  version  valid  to  70  Rg,  and  a  short  version  valid 
to  30  Rg. 

As  for  the  Mead-Fairfield  model,  parameter  sets  are  given  for  each  of  several  Kp  bins  (11  in  the 
Tsyganenko-Usmanov  model,  6  in  the  Tsyganenko  long  version,  and  8  in  the  Tsyganenko  short  version). 


2. 1.1. 5  Stern  Parabolic  Magnetopause  -  Stretched  Magnetosphere  Model 

The  field  due  to  the  magnetopause  currents  is  represented  within  this  boundary  as  the  gradient  of  a  scalar 
potential  satisfying  Laplace’s  equation  [Stern,  1985).  Ilie  requirement  of  a  closed  magnetosphere,  i.  e., 
that  the  field  lines  of  the  total  field  do  not  cross  this  boundary,  leads  to  the  condition  that  the  normal 
component  of  the  boundary  field  must  cancel  the  normal  component  of  the  total  field  due  to  internal 
sources  (dipole,  ring  current,  etc.).  Since  the  boundary  is  assumed  to  be  parabolic,  the  solution  is 
expanded  in  parabolic  harmonics  analogous  to  the  use  of  spherical  harmonics  for  problems  with 
spherical  boundaries.  The  axis  of  the  boundary  is  the  GSM  axis,  and  its  nose  is  at  local  noon.  No 
attempt  is  made  at  a  self-consistent  solution  of  the  pressure  balance  condition.  Thus,  there  is  no 
dependence  of  the  shape  on  the  dipole  tilt. 

The  ring  current  field  is  modeled  with  the  Tsyganenko-Usmanov  expression.  ITie  total  field  due  to  the 
dipole  and  cross-tail  currents  is  simulated  by  a  stretched  magnetosphere  model  [5fern,  1987).  The  field 
at  any  point  is  represented  in  terms  of  the  dipole  field  at  some  other  location,  specified  by  a  stretching 
function.  Both  the  stretched  dipole  and  ring  current  fields  are  accounted  for  in  the  boundary  conditions. 
The  solution  for  the  pure  dipole  has  a  simple  scaling  property  with  respect  to  pure  linear  scaling  of  the 
boundary.  The  ring  current  solution  has  a  similar  property,  if  the  dimension  of  the  ring  current  is  also 
scaled. 

2.1. 1.5.1  Model  Derivation 

The  basic  quiet  magnetopause  paraboloid  is  determined  by  the  requirements  that  the  noon  equatorial 
crossing  distance  (standoff  distance)  be  10  Rg,  and  the  dawn  and  dusk  equatorial  crossing  distances  be 
15  Rg  which  are  consistent  with  observations  [Fairfield,  1971  j.  ITie  expansion  for  the  tail  field  scalar 
potential  is  truncated  to  10  terms.  Ihe  dependence  of  the  magnetopause  field  on  standoff  distance  is 
derived  by  assuming  that  the  dimensions  of  the  magnetopause  vary  linearly  with  the  standoff  distance, 
and  then,  using  the  above-mentioned  line^ir  scaling  property  of  the  field. 
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For  the  ring  current,  there  is  no  specification  of  parameters.  We  have  assumed  strength  parameter 
BQ=-20nT  and  length  parameter  (Po)=4  R,;.  The  stretch  function  is  a  spline  fit  of  a  tabulation 
generated  by  the  author  and  is  designed  to  fit  night  side  equatorial  observations  of  Fairfield  [1968]. 


2.1.2  Computation  Times 

Table  3  lists  the  VAX  8650  computation  time  per  10000  subroutine  calls  for  the  various  external  models 
and  the  IGRF85  internal  model.  The  table  shows  that  the  time  for  a  total  internal  +  external  model 
calculation  is  not  terribly  sensitive  to  the  choice  of  external  model,  unless  the  external  model  is  Stern’s. 
Most  of  the  Stern  model  computation  time  is  in  the  magnetopause  part,  which  requires  generation  of 
Bessel  functions  for  the  parabolic  harmonic  expansions.  The  Bessel  functions  were  computed  with 
routines  from  Press,  et.  al.  [1986]. 


Table  3.  Vax  8650  Computation  Times  (sec/10000  calls) 


Olson-Pfitzer  Tilt  Dependent  2.0 

Olson-Pfitzer  Dynamic  5.2 

Mead-Fairfield  0.3 

Tsyganenko-Usmanov  1.8 

Tsyganenko  Short  Version  3.0 

Tsyganenko  Long  Version  3.6 

Stern  Parabolic  Magnetopause/Stretched  Magnetosphere  12.0 

IGRF85  9.2 


2.1.3  L  Parameter  Comparisons 

As  explained  by  Jordan  and  Bass  [1990],  the  Mcllwain  L  parameter  [McBwain,  1961],  at  a  point,  is  a 
function  of  the  first  two  adiabatic  invariants  of  particles  mirroring  at  that  point  provided  the  magnetic 
field  is  static.  For  internal  fields,  we  may  consider  this  to  be  true,  if  we  neglect  the  slow  secular  change 
known  to  be  taking  place.  In  Figures  33  and  34,  we  show  AL  =  L-Lj,  where  L^  is  the  dipole  L 
parameter,  r/sin^0,  with  0  the  magnetic  colatitude  and  r  the  radial  distance  of  the  point  in  R^. 

In  the  figures,  we  show  contours  of  constant  AL  for  the  IGRF85  model,  updated  to  1  January  1990.  The 
individual  contour  plots  are  at  designated  values  of  L^,  and  are  on  magnetic  latitude-local  time  grids. 
The  dipole  tilt  is  zero,  meaning  that  the  dipole  vector  is  normal  to  the  Earth-sun  line,  and  the  solar 
magnetic  longitude  (SLON)  is  zero.  The  contours  cover  only  locations  above  the  surface  of  the  Earth. 
Thus,  as  Lj  increases,  the  amount  of  latitude  space  covered  increases. 

Near  the  equator  AL  is  less  than  0.1  but  it  exceeds  this  value  in  the  South  Atlantic  Anomaly  and  at  high 
latitudes.  As  Lj  increases,  AL  =  1  contours  appear  at  the  high  latitudes,  which  were  below  the  surface 
at  lower  Lj.  One  remarkable  feature  that  appears  is  that  a  contour  is  nearly  stationary  beyond  the  Lj 
threshold  for  its  appearance;  thus.  AL  is  approximately  independent  of  L^. 
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Figure 


Figures  35  and  36  similarly  display  AL  for  a  field  composed  of  the  1GRF85  internal  model  and  the 
Olson-Pfitzer  tilt  dependent  external  model.  In  Figure  36,  we  have  marked  the  local  minima  and  maxima 
with  "L"  and  "H",  respectively,  indicating  their  values.  Now,  as  Lj  increases,  the  contour  pattern  begins 
to  deviate  substantially  from  the  IGRF85  baseline.  At  high  L^,  there  is  a  large  area  around  the  sub-solar 
point  with  AL  <  -1,  and  a  corresponding  area  around  the  anti-solar  point  with  AL  >  1. 

To  examine  this  effect  further,  and  in  particular,  to  more  easily  compare  the  models,  we  have  plotted, 
in  Figures  37  and  38,  A L  vs  L^  at  0°  and  20®  dipole  latitude,  respectively,  for  both  noon  and  midnight, 
and  for  0°  and  18°  dipole  tilt.  ITie  latitudes  and  radial  range  were  chosen  to  reflect  the  planned 
CRRES  orbital  parameter  ranges.  We  compare  in  these  figures  the  IGRF85  model  alone  ("NO  Ext”, 
solid  line)  with  results  obtained  for  IGRF85  -I-  Tsyganenko-Usmanov,  Tsyganenko  short  version  (’Tsyg” 
1987),  Mead- Fairfield,  and  Olson-Pfitzer  tilt-dependent  models.  The  constant  and  very  low  AL  for  the 
internal  model  alone  is  in  agreement  with  what  was  observed  in  the  contours  plots.  One  readily  sees  the 
growing  effect  of  the  external  models  as  L^j  increases.  A  radiation  belt  particle  which  mirrors  at  the 
equator  (equatorial  pitch  angle  =  90°)  must  follow  a  path  of  constant  L.  At  noon  the  L=Lj  path  is 
outside  the  L^  =  constant  circle  by  approximately  the  amount  shown  in  the  noon  plot.  Thus,  at  noon, 
the  L=7  path  is  approximately  0.5  Rj;  outside  the  Lj=7  circle.  At  midnight,  nearly  the  opposite  is  true: 
the  L=7  path  is  inside  the  Lj=7  circle.  Ibus.  a  trapped  particle,  conserving  its  adiabatic  invariants  and 
mirroring  at  the  equator,  travels  a  path  significantly  distorted  from  a  circular  path. 

The  differences  beh^’een  the  models  is  up  to  0.25  Rj;  in  the  high  L^  region  of  the  plots.  The 
Olson-Pfitzer  tilt-dependent  model  yields  the  largest  L  values  in  the  inner  region,  while  the 
Mead-Fairfield  model  gives  the  largest  values  in  the  outer  region  at  noon,  and  the  Tsyganenko  and 
Tsyganenko-Usmanov  models  yield  the  largest  values  at  midnight.  The  Olson-Pfitzer  midnight  L 
increases  nearly  linearly  with  distance,  while  the  L  values  for  the  other  models  increase  more  rapidly  in 
the  outer  region  than  in  the  inner  region.  The  dependence  of  the  L  value  on  dipole  tilt  is  negligible  at 
the  equator,  but  significant  at  20°  latitude.  Thus,  we  see  that  there  are  substantial  modifications  in  the 
trapped  particle  orbits  when  the  external  field  is  introduced  and  noticeable  differences  between  various 
models  of  the  external  field.  Thus,  it  is  worthwhile  to  compare  the  models  with  observations. 


2.1.4  Models  vs.  Data 

A  static  radiation  belt  model  will  require  a  good  static  magnetic  field  model.  ITierefore,  we  have  chosen 
SCATHA  magnetometer  data  from  the  magnetically  quiet  20  April  1979  for  comparison  with  the  models. 
Despite  the  quiet  condition,  the  solar  wind  pressure  increased  steadily  during  the  period,  allowing  us  to 
test  the  magnetopause  modeling  capabilities  of  the  various  models.  ITie  following  day,  there  was  a  small 
substorm,  during  which,  the  measured  magnetic  field  magnitude  decreased  considerably  from  the 
prediction  of  the  static  Olson-Pfitzer  tilt  dependent  model.  If  we  wish  to  model  adiabatic  variations  of 
the  radiation  belt  particles  during  such  substorms,  we  need  to  have  an  accurate  model  of  the  magnetic 
field  variations. 

Since  the  SCA I  HA  orbit  is  limited  to  the  outer  belt,  we  also  examined  some  published  AMP'lE/CCE 
results  [Fairfield,  et.  al.,  1987).  As  argued  by  these  authors,  the  inner  magnetic  flux  has  a  bearing  on 
the  equatorial  crossings  of  magnetic  field  lines  intersecting  the  Earth’s  surface  at  given  latitudes.  Ilie 
flux  intersecting  the  surface  equatorward  of  a  field  line  in  a  specified  longitudinal  sector  is  approximately 
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equal  to  the  flux  intersecting  the  equator  in  the  same  sector  Earthward  of  the  equatorial  intersection  of 
the  field  line. 


2. 1.4.1  Comparisons  of  Models  with  SCATHA  Data 

The  SCATHA  satellite  orbit  is  a  5.3  Rg  x  7.8  Rg  low  inclination  (7.9°)  near-geosynchronous  orbit  with 
an  easterly  drift  rate  of  about  5°  per  day.  On  20  April  1979,  the  perigee  of  the  SCATHA  orbit  was  at 
1600  hours  MLT,  the  apogee  was  at  0300  hours  MLT,  and  the  magnetic  latitude  was  within  ±4.5°.  The 
magnetometer  data  used  in  the  following  were  averaged  over  1/20  Rg  bins  (-15  min.). 

Figure  39  shows  the  differences  between  the  measured  and  modeled  fields  for  this  day.  The  internal  field 
model  used  is  the  Barraclough  1975.  updated  to  the  date  in  question. 

Of  all  the  models,  the  Stern  model  disagrees  the  most  with  the  data,  while  the  Olson-Pfitzer  dynamic 
model  agreement  is  generally  the  best,  particularly  at  perigee,  where  the  differences  are  greatest.  With 
the  exception  of  the  Stern  and  Olson-Pfitzer  dynamic  models,  the  differences  at  the  end  of  the  day,  as 
the  next  perigee  is  approached,  are  positive  but  opposite  in  sign  from  the  differences  at  the  beginning 
of  the  day  at  the  same  location  in  the  orbit.  Thus,  the  field  seems  to  have  increased  over  the  24  hour 
period  relative  to  the  models.  In  fact,  the  standoff  distance  implied  by  the  solar  wind  data  decreased 
steadily  during  the  day  (Figure  40).  Thus,  it  appears  that  the  Olson-dynamic  and  the  Stern  models,  with 
their  explicit  dependence  on  the  standoff  distance,  depict  this  variation  better  than  do  the  other  models. 

The  comparatively  large  deviations  of  the  Stern  model  are  not  surprising  since  that  model  has  not  been 
extensively  fit  to  data.  The  large  negative  deviations  for  the  day  local  times  (meaning  the  model  field 
intensity  exceeds  the  measured)  could  possibly  be  improved  by  increasing  the  ring  current  but  this  would 
worsen  the  disagreement  of  opposite  sign  at  night  local  times.  The  stretch  function  evidently 
overestimates  field  line  stretching  for  this  case  of  extremely  low  activity. 

Figure  41  shows  the  measured  model  magnetic  field  differences  for  the  next  day,  21  April  1979.  On  that 
day,  a  substorm  occurred  with  growth  phase  commencing  at  -  5  hours  U 1.  The  measured  magnetic  field 
is  seen  to  decrease  significantly  with  respect  to  all  the  models.  Initially  the  Olson-Pfitzer  dynamic  model 
seems  to  follow  the  data  marginally  better  than  the  other  models,  but  that  may  only  be  because  it  uses 
hourly  indices,  while  the  other  dynamic  models  use  the  3-hourly  Kp.  TTius.  later  on,  the  others  are  doing 
just  as  well.  The  start  of  the  substorm  is  indicated  by  the  sharp  drop  in  the  measured  minus 
Olson-Pfitzer  static  model  difference  (solid  line)  at  -  4-5  hours  UT.  During  this  period,  the  measured 
minus  Olson-Pfitzer  dynamic  model  difference  (dotted  line)  actually  increases,  meaning  that  the  dynamic 
model  depicted  a  faster  decrease  than  was  observed.  Later,  however,  the  Olson-Pfitzer  dynamic  model 
is  almost  as  far  off  the  data  as  the  others.  During  subsequent  disturbances  on  that  day,  no  model  does 
outstandingly  well. 
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Figure  40.  Noon  magnetopause  standoff  distance:  20  April  1979. 
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2. 1.4.2  Midnight  Equatorial  Field  Strength  Depression 


One  important  feature  of  the  magnetic  field  is  the  variation  of  the  field  magnitude  along  the  equator  in 
the  noon-midnight  plane.  On  the  midnight  side,  the  external  sources  generally  cause  a  depression  in  the 
total  field  strength,  stretching  tailward  the  field  lines  connected  to  given  locations  on  the  Earth’s  surface, 
llierefore,  as  pointed  out  by  Fairfield,  et.  at.  [19871,  the  magnetic  mapping  of  given  points  on  the 
Earth’s  surface  to  the  equator  is  determined  by  the  amount  of  this  depression.  Thus,  an  accurate  model 
of  the  midnight  equatorial  depression  is  needed  to  relate  auroral  observations,  such  as  precipitation 
boundaries  and  field  aligned  currents,  to  high  latitude  equatorial  factors,  such  as  convection,  injections, 
and  plasma  pressure  gradients. 

Figure  42  compares  averaged  observed  equatorial  midnight  field  depressions,  adopted  from  Fairfield,  et. 
al.  [1987],  for  various  Kp  bins  with  those  predicted  by  the  models.  The  0-1  bin  and  the  2-3  bin  are  from 
Sugiura  and  Poros  [1973],  based  on  OGO  measurements,  while  the  results  for  the  other  two  bins  were 
determined  by  Fairfield,  et.  al.  [1987]  from  the  AMPTE/CCE  data.  According  to  the  latter  authors, 
some  differences  in  the  processing  of  the  data  may  have  led  to  differences  in  the  high  Kp  results  for  the 
two  satellites  at  the  larger  distances;  however,  the  smaller  distance  values  seem  consistent. 

For  each  model  the  quantity  plotted  is 

♦  B„  I  -  Pp.p  1  <“l 


for  zero  dipole  tilt,  'fhe  Olson-Pfitzer  tilt-dependent  model  is  presented  only  for  the  Kp  =  0-1  bin.  For 
the  Mead-Fairfield,  Tsyganenko-Usmanov,  and  Tsyganenko  1987  short  models,  we  chose  bins  as  close 
to  those  of  the  data  as  the  binning  structure  of  the  individual  models  permitted.  Thus,  the 
Nfead-Fairfield  model  is  specified  for  4  bins.  Of  these,  only  two.  Kp  <  2  and  K  >  3  were  sufficiently 
comparable  to  the  data  bins.  For  the  Tsyganenko-Usmanov  and  Tsyganenko  Iw  models,  the  results 
were  obtained  by  weighted  averaging  of  appropriately  selected  model  bins;  the  number  of  obser\ations 
in  each  such  bin,  given  by  the  authors,  determined  the  weighting  factor  for  that  bin. 

For  the  Olson-Pfitzer  model,  the  mean  standoff  distance  and  value  for  each  K  bin  w'ere  determined 
from  1983  data  on  the  NSSDC  OMNI  online  data  base.  ITiese  are  given  in  Table  4.  In  the  Stern  model, 
as  for  the  Olson-Pfitzer  tilt  dependent  model,  only  a  single  set  of  results  is  given. 


Table  4.  Average  Values  for  the  Olson-Pfitzer  Dynamic  Model 


Standoff  Distance 

Ds, 

No.  of  points 

0-1 

9.680 

-1.1 

440 

0-2- 

9.601 

-2.3 

819 

2-3 

9.091 

-7.9 

1098 

>  3  + 

8.514 

-28.0 

1262 
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Figure  42.  Averaged  measured  [Fairfield,  et.  al.,  1987]  and  model  midnight  equatorial  magnetic  field  depression. 

'Hie  K  -  0-1  and  2-3  measurements  are  from  OGO  data,  the  K  <  2  and  K  >  3+  measurements  are  from  AMPTE  data. 


ITie  Olson-Pfitzer  tilt-dependent  model  agrees  quite  well  with  the  measurements  for  the  low  Kp  bin. 
This  is  not  surprising,  as  this  model  is  also  derived  from  OGO  data.  The  Mead-Fairfield, 
Tsyganenko-Usmanov.  and  Tsyganenko  models  all  predict  depressions  smaller  than  the  data  at  smaller 
distances.  The  Kp  dependence  in  the  Mead-Fairfield  model  is  nearly  constant  with  distance,  while  in  the 
other  two  models  the  Kp  dependence  decreases  with  increasing  distance,  in  apparent  contrast  with  the 
data.  The  disagreement  at  small  distances  is  probably  due  to  the  lack  of  equatorial  data  at  these 
distances  in  the  data  sets  used  to  derive  these  models,  llie  Olson-  Pfitzer  dynamic  model  exhibits 
considerably  greater  depression  than  the  data.  Like  the  Isyganenko-Usmanov  and  Tsyganenko  model 
depressions,  that  for  the  Olson-Pfitzer  dynamic  model  exhibits  decreasing  Kp  dependence  with  increasing 
distance.  The  Stern  model  seems  to  fit  the  Kp  <  2-  measured  depression,  which  is  not  too  surprising, 
as  Stern’s  dipole  stretch  model  is  based  on  data  averaged  over  all  geomagnetic  activity  levels. 


2.1.5  Summary 


Six  models  of  the  external  magnetospheric  field  have  been  reviewed.  The  purpose  of  these  models  is  to 
provide  accurate  and  computationally  efficient  evaluation  of  the  contribution  to  the  magnetospheric 
magnetic  field  due  to  external  sources,  for  use  in  analysis  and  interpretation  of  magnetospheric  data. 

Their  properties  were  examined  in  the  .5-8  Rj;  range,  that  portion  of  the  CRRES  orbit  where  the  external 
field  is  expected  to  be  important.  An  external  field  typically  modifies  the  L  parameter  beyond  4  Rg. 
The  dependence  of  L  on  the  dipole  tilt  is  small  along  the  dipole  equator,  but  significant  at  higher 
latitudes. 

For  a  magnetically  quiet  period,  the  magnetic  field  magnitudes  derived  from  all  the  models  agree  closely 
with  that  derived  from  SCATHA  magnetometer  observations  over  that  period.  Upon  closer  review,  it 
was  found  that  the  Olson-Pfitzer  dynamic  model  tracks  the  data  the  best,  since  its  magnetopause  model 
depends  explicitly  on  the  standoff  distance.  Stern’s  model,  although  its  magnetopause  model  is  also  tied 
to  the  standoff  distance,  showed  the  poorest  agreement  with  the  data,  possibly  because  its  magnetopause 
model  is  purely  theoretical,  and  its  tail  field,  implied  by  the  stretch  model,  is  too  large  (and  negative)  on 
the  night  side  and  identiailly  zero  on  the  day  side.  During  a  substorm  on  the  following  day,  none  of  the 
models  predict  the  large  drop  in  the  field  magnitude  that  is  seen  by  the  data,  although  several  of  the 
models  show  a  very  small  decrease. 

For  zero  tilt,  the  Stern  and  Olson-Pfitzer  tilt-dependent  model-calculated  equatorial  midnight  field 
depressions  are  consistent  with  averaged  AMPTE  and  OGO  measurements  as  functions  of  distance.  The 
Olson-Pfitzer  dynamic  model  predicts  larger  depressions  than  observed,  while  the  other  dynamic  models 
predict  smaller  depressions. 
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2.2.1 


After  an  initial  study  of  the  available  external  magnetospheric  magnetic  field  models,  we  found  two  new 
models  which  had  not  been  covered:  Tsyganenko  1989  and  Hilmer- Voigt.  Thus,  a  subsequent  study  was 
performed  to  better  evaluate  these  models  along  with  the  Olson-Pfitzer  Dynamic  (1988)  model.  In  fact, 
two  similar  studies  were  performed,  the  first  was  a  comparison  to  SCATHA  data  (20  and  21  April  1979) 
and  the  second,  a  comparison  to  CRRES  data  (26  August  1990).  Tsyganenko  1989  started  with  an 
axisymmetric  infinitely  thin  current  disk  for  the  tail  current.  He  transformed  the  disk  to  allow  for 
hinging  (where  the  tail  current  bends  away  from  the  dipole  equatorial  plane  in  the  near  Earth  region 
to  become  parallel  to  the  plasma  sheet  in  the  GSM  x-y  plane).  Also,  the  disk  was  modified  to  allow  for 
transverse  bending  in  the  y-direction.  Factors  were  added  to  drop  off  smoothly  as  the  magnetopause 
was  approached.  The  ring  current  is  represented  by  the  second  derivative  of  the  axisymmetric  infinitely 
thin  current  disk  to  better  confine  its  contribution  to  the  field.  Only  westward  current  was  represented. 
There  is  no  magnetopause.  There  is  no  inner  edge  to  the  current  sheet.  No  data  was  included  inside 
of  four  Rg.  Thus,  this  model  is  not  really  intended  for  use  in  the  inner  magnetosphere.  Finally,  the  data 
was  binn^  by  Kp  and  fit  accordingly.  Ihus,  this  is  an  average  model  of  the  field  within  each  Kp  bin. 
As  such,  we  found  it  to  be  in  relatively  good  agreement  with  the  data.  It  is  not  suited  for  looking  at 
dynamic  changes  in  the  data,  such  as  the  compression  of  the  field  as  the  solar  wind  pressure  increases 
and  the  magnetopause  moves  in.  Nor  can  it  resolve  fine  scale  activity  (Kp  is  only  available  as  three  hour 
averages).  However,  as  an  average  model,  it  did  reasonably  well  in  comparison  to  both  our  quiet  time 
and  early  storm  time  data  (Figures  4.'!  and  44).  Ifs  fast  and  very  easy  to  use. 

2.2.2  Hilmer-Voiet  Model 

Hilmer-Voigt  is  driven  by  four  physical  input  parameters:  dipole  tilt,  stand-off  distance  of  the 
magnetopause,  activity  index  Dst,  and  the  midnight  equatorward  boundary  of  the  diffuse  aurora. 
Considering  all  of  the  models  that  we  have  looked  at,  this  is  the  most  flexible.  The  magnetopause  is 
represented  by  Voigt’s  1981  shielded  vacuum  dipole  model  with  a  fixed  shape  of  a  hemisphere  attached 
to  a  cylinder.  The  ring  current  is  represented  by  a  vector  potential  similar  to  that  used  by  Tsyganenko 
1987,  but  Hilmer’s  expression  includes  eastward  travelling  current.  The  tail  current  is  based  on  a  method 
used  by  Tsyganenko  and  Usmanov  which  uses  magnetic  filaments  extending  infinitely  in  y.  The  filaments 
are  grouped  in  segments.  Transverse  bending  and  hinging  are  included.  The  stand-off  distance 
determines  the  strength  of  the  magnetopause  currents.  Dst  controls  the  ring  current,  and  the  midnight 
equatorward  boundary  maps  to  the  inner  edge  of  the  tail  current  sheet. 

Therefore,  this  model  is  the  most  flexible  in  rendering  the  dynamic  changes  of  the  magnetosphere. 
However,  it  is  also  the  slowest  and  most  difficult  to  use.  The  four  input  parameters  arc  used  to  find 
intermediate  parameters  which  are  actually  required  to  run  the  model.  Until  very  recently,  this  input 
parameter  selection  code  was  not  available,  so  one  had  to  contact  Hilmer  to  get  the  necessary  parameters 
to  run  his  code  for  a  given  configuration  of  the  field.  We  now  have  a  preliminary  release  of  this 
intermediate  code  to  work  with  and  test.  Nonetheless,  this  model  agreed  best  with  the  data  overall 
(Figures  43  and  44).  It  is  not  really  suitable  for  heavy  data  processing,  but  is  very  useful  for  modeling 
dynamic  change  on  a  fairly  small  saile. 
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B-ne!d  MagnUuds  CRRES  Orbit  76  26  August  1990 
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Figure  43.  lliis  is  a  sample  of  a  preliminary  comparison  of  these  models  to  CRRES  data.  'Hie  top 
panels  show  the  data  and  model  |.H|  for  part  of  26  August  1990.  ITie  bottom  panel  shows  the  difference 
ratios  seen  between  the  models  and  data  for  \B\. 
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3.1  Magnetic  Latitude  and  Local  Time  at  Bmin 


Along  a  field  line  in  the  trapped  radiation  belts,  there  exists  a  point  of  minimum  magnetic  field  strength. 
All  particles  bouncing  along  the  field  line  must  pass  though  this  minimum.  Thus,  all  trapped  particles 
in  the  radiation  belt  must  pass  through  this  surface  of  minima  which  we  often  refer  to  as  BMIN. 
Therefore,  from  a  model  of  the  trapped  populations  along  this  surface,  it  is  possible,  in  principle,  to 
derive  the  populations  at  all  other  points.  In  practice,  however,  such  a  derivation  requires  the 
determination  of  the  BMIN  location  and  properties  as  a  function  of  time  and  position.  This  requires 
that  we  find  the  intersection  of  the  field  line  through  the  starting  position  with  the  BMIN  surface.  Since 
straightforward  field  line  tracing  can  be  quite  expensive  computationally,  we  have  done  a  survey  of  BMIN 
surface  properties  in  the  hope  of  finding  properties  that  would  lead  to  a  simpler  procedure. 

Ideally,  the  BMIN  surface  should  coincide  with  the  dipole  equator.  However,  it  is  known  that  at  large 
distances  in  the  midnight  local  time  sector,  the  solar  wind  alters  the  surface  significantly  when  the  tilt 
of  the  dipole  with  respect  to  the  earth-sun  line  differs  from  zero.  In  such  cases,  the  dipole  equator  is 
not  parallel  to  the  solar  wind.  At  small  distances  from  the  earth,  the  internal  field  dominates.  Thus,  the 
minimum  still  lies  very  close  to  the  equator.  As  the  distance  increases,  the  solar  wind  forces  the  BMIN 
surface  to  curve  away  from  the  dipole  equator,  toward  a  surface  parallel  to  the  earth-sun  line.  In 
addition,  it  causes  magnetic  field  lines  to  be  significantly  twisted  in  local  time,  so  that  lines  emanating 
from  the  earth’s  surface  at  dawn  or  dusk  local  time  are  twisted  toward  midnight.  Thus,  the  local  time 
at  a  given  CRRES  location  may  differ  from  the  local  time  at  the  appropriate  BMIN  intersection. 

To  study  these  effects  with  CRRES,  we  have  constructed,  from  the  standard  ephemeris  files  for  the  first 
450  orbits,  a  BMIN  survey  data  base  which  contains  for  each  ephemeris  point: 

orbit  number  and  time  (year,  day  number,  UT) 

solar  magnetic  coordinates  (radial  distance,  latitude,  local  time)  at  CRRES 
dipole  tilt 

magnetic  longitude  of  the  sun 

solar  magnetic  coordinates  at  BMIN 

magnetic  field  strength  and  L  at  both  BMIN  and  CRRES 

Figure  45  shows  BMIN  local  time  and  orbit  number  for  all  BMIN  distances  greater  than  7  R^^.  It  shows 
that  apogee  approaches  local  midnight  just  prior  to  the  end  of  this  period.  ITie  maximum  solar  magnetic 
latitude  magnitude  on  the  BMIN  surface  was  4.2°.  The  maximum  local  CRRES-BMIN  local  time 
difference  was  11  minutes.  Figure  46  shows  the  BMIN  distance  and  solar  magnetic  latitude  for  BMIN 
distances  greater  than  6  Rj;.  Although  there  is  a  lot  of  variation  in  the  data,  there  is  a  trend  toward 
increasing  northern  latitude  with  distance.  The  apogee  approaches  midnight  in  midwinter,  which  is 
dominated  by  high  negative  dipole  tilt.  In  this  situation  the  BMIN  surface  curves  northward  of  the 
dipole  equator. 

In  conclusion,  the  solar  wind  does  not  distort  the  minimum  B  surface  severely  from  that  expected  for 
the  simple  dipole  model,  nor  does  the  local  time  vary  significantly  along  a  field  line.  Such  variations  that 
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Figure  45.  Solar  magnetic  (SM)  local  time  and  orbit  number  for  all  cases  of  BMIN  distance  greater  than 
7  Rg,  orbits  1  -  450. 
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Figure  46.  SM  latitude  and  radial  distance  for  all  cases  of  BMIN  distance  greater  than  6  R^.  orbits  1  - 
450. 
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do  exist  should  be  modeled  analytically  to  estimate  the  intersection  of  a  field  line  with  the  BMIN  surface., 
as  an  alternative  to  direct  tracing. 


3.2  B/Bp  vs  L  and  Magnetic  Latitude 


Existing  software  has  been  modified  and  improved  to  compute  the  relative  dwell  times  of  a  satellite  for 
a  specified  range  of  B/B^  and  L-shell  values  over  a  specified  time  interval.  The  dwell  chart  provides  a 
graphical  representation  of  the  fraction  (percentage)  of  time  that  a  satellite  will  spend  in  a  given  B/B^ 
and  L-shell  interval.  The  data  is  represented  either  by  a  black  and  white  density,  or  by  different  colors. 
A  color  Dwell  Chart  was  produced  for  the  CRRES  satellite  for  a  one  year  period  using  updated  CRRES 
orbital  elements. 

A  color  Dwell  Chart  was  produced  for  the  CRRES  satellite  for  a  1/2  year  period  (July  1990  -  January 
1991)  for  the  researcher  to  supplement  the  1  year  color  Dwell  chart  produced  earlier.  Software  was 
written  to  produce  Black  and  White  versions  of  the  same  Dwell  Charts  using  the  HPGL  graphics 
capabilities  of  the  HP  LaserJet  III  printer.  Figure  47  shows  the  black  and  white  version  of  the  1  year 
Dwell  Chart. 

Code  was  developed  which  would  read  the  CRRES  orbit  database  files  and  compute  the  time  of  L  value 
crossings,  the  value  of  B/B^,  and  the  magnetic  latitude  (SM.  GM  and  GSM  coordinates)  at  the  crossing 
times.  A  series  of  scatter  plots  (B/B^,  vs.  magnetic  latitude  (SM,  GM  and  GSM))  for  selected  L  values, 
were  produced  for  orbits  1  -  510.  Figure  48  is  an  example. 
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4.0  INSTRUMENT  MODELING 


4.1  PROTEL  Contamination  Q>de 


4.1.1  Motivation 

The  Proton  Telescope  (PROTEL)  is  one  of  the  instruments  on  the  CRRES  satellite.  A  detailed 
description  of  the  PROTEL  instrument  is  provided  in  a  recent  report  [Lynch,  et  aL,  1989].  PROTEL 
measures  proton  counts  in  the  1  to  100  MeV  energy  range  in  the  radiation  belts.  It  consists  of  two 
detector  instruments,  the  low  energy  head  [LEH]  (1-9  MeV)  and  the  high  energy  head  [HEH]  (6  -  100 
MeV),  and  a  dedicated  processor  which  processes  the  raw  data  from  the  detectors,  and  at  1  second 
intervals,  transfers  the  reduced  data  to  the  satellite’s  telemetry  system.  The  reduced  data  consists  of 
counts  in  24  energy  channels  (8  for  LEH,  16  for  HEH)  spaced  logarithmically  in  the  1  - 100  MeV  energy 
interval,  together  with  environmental  data  and  raw  counts  for  the  solid  state  particle  detectors. 

PROTEL  is  designed  to  operate  in  a  hostile  environment,  where  it  is  subjected  to  a  high  density  flux  of 
electrons,  protons  and,  to  a  lesser  extent,  heavier  ions  and  cosmic  radiation.  The  design  integrates 
passive  shielding  techniques,  magnetic  deflection  of  electrons,  and  detection  (coincidence  and 
anticoincidence)  logic  which  is  used  to  eliminate  spurious  proton  counts. 

In  order  to  properly  evaluate  the  performance  of  PRO  lEL,  it  is  necessary  to  model  the  behavior  of  the 
PROIEL  in  its  hostile  environment.  In  order  to  do  so,  it  is  necessary  to  anticipate  deviations  of  the 
PROTEL  instrument  from  its  nominal  design,  and  to  identify  the  potential  sources  of  errors  due  to 
limitations  of  the  PRO'lEL  design  and  due  to  the  hostile  environment.  In  the  context  of  this  report, 
contamination  refers  to  the  penetration  of  high  energy  protons  from  angles  outside  the  nominal 
acceptance  cone  which  falsely  trigger  the  PROIEL  detection  logic,  despite  the  use  of  passive  shielding, 
consisting  of  aluminum,  brass,  iron  and  tungsten  in  the  construction  of  the  PROTEL  housing. 

In  a  recent  report  [Redus,  et  al.,  1990],  it  has  been  suggested  that  PROTEL  will  have  problems 
resulting  from  the  penetration  of  high  energy  protons  from  certain  angles  outside  the  entrance  cone. 
These  protons  will  produce  false  counts  (contamination)  in  the  data.  The  effect  of  contamination  on 
observed  counts  will  strongly  depend  upon  the  energy  spectrum  and  angular  distribution  of  the  high 
energy  protons.  For  this  reason,  a  software  package  (PRO T  EL  Contamination  Code)  which  models  the 
contamination  problem  for  the  HEH  has  been  developed. 

Portions  of  the  PROTEL  contamination  code  have  been  used  to  model  the  PROTEL  calibration 
experiments  which  were  performed  using  an  accelerator  at  Harvard  University.  The  predictions  of  the 
contamination  code  were  compared  with  the  Harvard  experimental  data.  Within  the  limits  of 
experimental  error,  the  agreement  between  the  predictions  and  the  model  appears  to  be  satisfactory 
[Hein.  1990). 
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4.1.2  Methodology 


The  PROTEL  HEH  contains  six  solid  state  detectors  D1-D6.  D1-D5  each  contain  two  detection  areas; 
the  first  is  a  disk,  and  the  second  is  a  ring  surrounding  the  disk.  D6  has  only  a  circular  disk.  Pulse 
height  circuitry,  coincidence  and  anti-coincidence  circuitry,  and  logic  are  used  to  classify  detection  events 
in  the  several  particle  detectors.  A  detection  event  is  recorded  in  one  of  the  energy  channels  only  if  it 
meets  the  energy  deposit  requirements  of  the  detectors  associated  with  the  energy  channel  and  the 
requirements  of  the  coincidence/anticoincidence  logic.  The  ring  detectors  and  D6  are  used  in  connection 
with  the  anti-coincidence  logic 

4.1.2.1  Physical  Assumptions 

The  physical  assumptions  made  in  developing  the  PROTEL  contamination  code  are  as  follows: 

The  primary  effect  is  due  to  protons  which  penetrate  the  PROIEL  housing/shielding  and  are  decelerated 
in  the  direction  of  motion,  while  interacting  with  the  electrons  in  the  material  through  which  they  pass. 
The  protons  are  assumed  to  travel  in  straight  lines  until  they  either  are  stopped  by  or  penetrate  the 
material  they  are  passing  through.  This  is  a  relatively  good  assumption  for  high  energy  protons. 
However,  important  effects  are  neglected,  such  as  Coulomb  scattering  with  the  nuclei,  inelastic  scattering 
off  of  the  nuclei,  including  “SI  ARS"  -  nuclear  reactions  which  produce  protons,  neutrons  and/or  gamma 
rays. 

These  assumptions  permit  the  modeling  of  contamination  by  means  of  "ray  tracing"  through  the 
PROTEL  HEH  for  particles  entering  from  a  given  direction  and  with  a  given  energy  as  described  by  a 
differential  flux.  The  Harvard  Accelerator  experiments  were  modeled  assuming  that  the  accelerator  beam 
was  a  parallel  beam  with  little  variation  transverse  to  the  direction  of  the  beam.  In  the  more  general 
case,  it  was  assumed  that  the  energy  spectrum  and  the  angular  distribution  of  the  incident  radiation  are 
independent.  A  collection  of  representative  points  on  the  first  detector  is  selected,  together  with  a 
direction.  For  each  point,  the  spectrum  is  computed  (using  the  Janni  energy/range  relation)  at  the  point 
of  the  first  detector.  A  determination  is  then  made  as  to  which  detectors  the  ray  passes  through;  the 
corresponding  energy  deposited  in  each  detector  is  calculated,  and  the  corresponding  channel  counts  are 
computed,  llie  counts  over  the  detector  surface  for  a  given  direction  are  summed,  and  the  sum  is 
weighted  according  to  the  direction.  The  weighted  channel  counts  are  summed  to  obtain  a  total  count 
for  the  input  spectra. 

4. 1.2.2  Janni  Ranee-Energy  Tables 

The  Janni  Range-Energy  Tables  jJanni,  1982)  were  incorporated  into  the  PROTEL  Contamination  Code 
for  the  following  materials:  Aluminum,  Brass,  Iron.  Tungsten  (for  the  PRO'lliL  HEH  passive  shielding 
and  housing),  Silicon  (for  the  particle  detectors)  and  Lucite  (for  the  absorbers  used  in  the  Harvard 
accelerator  experiments). 

4. 1.2.3  Rav  Tracing  of  PROTEL  HEH 

FORTRAN  Code  was  developed,  which  performed  a  ray  trace  through  the  different  materials  present 
in  the  PROTEL  HEH  passive  shielding  and  housing.  Careful  measurements  were  made  of  the  original 
blueprints  of  the  different  materials  used  in  the  PR01E.L  HEH  housing  and  passive  shielding.  A 
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routine  was  developed,  which  computes  the  thicknesses  of  the  different  materials  through  which  a  proton 
would  pass  as  it  entered  from  a  given  direction.  Ihe  routine  would  return  arrays  which  identified  the 
type  of  material,  and  the  thickness  encountered  to  an  accuracy  of  0.1  mm.  The  latter  was  used  to 
compute  the  proton’s  energy  loss  using  the  Janni  Range-Energy  Relation. 

4.1.2.4  Monte  Carlo  Technique 

Early  in  the  development  of  the  Protel  Contamination  Code,  it  was  discovered  that  the  computational 
overhead  of  doing  a  standard  (deterministic)  computation  was  prohibitive.  For  the  Harvard  experimental 
data  computation,  only  a  single  direction  was  required,  and  the  computation  effort  required  was 
reasonable.  For  the  full  contamination  computation,  it  was  desired  to  use  Monte  Carlo  methods,  in 
effect,  to  perform  an  integration.  For  each  angular  distribution  of  interest  (isotropic,  and  two  cases  of 
mirror  plane  geometry  corresponding  to  maximum  and  minimum  obstruction  due  to  the  electron 
deflection  magnet,  and  the  sin^  a  pitch  angle  distribution)  code  was  written  in  which  20,000  points  on 
the  first  detector,  and  a  direction  (appropriate  to  the  angular  distribution)  were  randomly  selected.  The 
tracing  was  then  performed,  which  provided  a  table  of  materials  and  thicknesses  encountered  for  each 
point.  The  response  function  was  then  computed  for  proton  energies  ranging  from  5-110  MeV  at  0.1 
MeV  intervals,  and  from  1 10-250  MeV  at  1  MeV  intervals.  Ihe  response  function  could  then  be  used 
to  compute  the  results  for  each  angular  distribution  for  an  arbitrary  energy  distribution.  In  order  to 
verify  the  sufficiency  of  20.000  points,  comparison  runs  were  made  using  different  "seeds"  for  the  random 
number  generation,  with  the  result  that  the  output  agrees  within  a  few  percentage  points. 


4.1.3  Results 

Ihe  response  function  for  Channel  10  is  provided  for  the  isotropic  case  (Figure  49)  and  for  the  Mirror 
Plane  Distribution  Case  (Figure  50).  For  the  mirror  plane  case,  two  curves  are  provided  which 
correspond  to  orientations  of  the  mirror  plane  relative  to  the  maximum  (0°)  and  minimum  (90°) 
obstruction  due  to  the  electron  deflection  magnet.  Note  that  the  effects  of  contamination  for  the  mirror 
plane  ease  are  not  as  severe  as  that  found  for  the  isotropic  case. 

4.1. 3.1  Harvard  Accelerator  Experiment 

Early  in  the  development  of  the  PRO  TEL  HEH  Contamination  Code,  portions  of  the  code  were  used 
to  compute  the  contamination  associated  with  the  PROTEL  calibration  tests  performed  with  the  Harvard 
cyclotron.  The  PROTEL  instruments  were  placed  at  several  different  orientations  at  several  different 
energies,  and  subject  to  proton  beams  of  different  energies,  where  the  latter  were  obtained  by  using 
different  thicknesses  of  lucite  as  a  beam  degrader.  Ihe  PROTEL  HEH  mass  distribution  model  and  the 
Janni  range  -  energy  relation  were  used  to  compute  the  counts,  assuming  that  the  proton  beam  was 
parallel  and  uniform  across  the  HEH.  Within  experimental  errors  of  the  calibration  experiment  and 
the  positioning  of  the  HEH,  the  agreement  was  reasonable. 
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Figure  49,  PROTEL  HEH  Response  Function  for  the  isotropic  case  for  channel  10.  The 
nominal  response  occurs  near  30  MeV.  TTie  contamination  is  significant  above  100  MeV,  and 
corrections  for  contamination  need  to  be  made. 
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Figure  50.  Response  Function  for  Mirror  Plane  Case,  for  two  orientations  of  the  mirror  plane 
with  respect  to  the  PROTEL  electron  deflection  magnet. 
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4.1.3.2  Contamination  Counts 


Predictions  from  the  Protel  Contamination  Code  have  been  produced  in  tabular  form,  and  also  in 
graphical  form  and  have  been  made  available  to  researchers  for  use  in  the  evaluation  of  the  performance 
of  the  PROTEL  HEH.  In  Table  5  the  response  function  was  used  to  compute  the  fraction  of  counts 
due  to  contamination  for  the  PRO  lEL  HEH  head  for  power  law  spectra  and  an  isotropic  distribution, 
llie  results  indicate  that  the  contamination  is  severe  for  flat  and  relatively  flat  spectra. 

Table  5 


FRACTION  OF  PROTEL  HEH  COUNl^S  DUE  TO  CONTAMINATION 

Isotropic  Distribution 
POWER  LAW  SPECTRUM  E*’ 


CHAN 

o 

II 

O' 

-1 

-2 

-3 

-4 

-5 

-6 

-7 

-8 

-9 

-10 

1 

.9889 

.7414 

.1358 

.0482 

.0355 

.0289 

.0241 

.0203 

.0172 

.0147 

.0125 

2 

.9751 

.6176 

.1057 

.0385 

.0263 

.0198 

.0154 

.0121 

.0095 

.0076 

.0060 

3 

.9534 

.5284 

.1672 

.1174 

.1035 

.0949 

.0882 

.0826 

.0778 

.0736 

.0697 

4 

.9513 

.5189 

.0982 

.0365 

.0231 

.0164 

.0122 

.0093 

.0071 

.0056 

.0044 

5 

.9686 

.6816 

.2003 

.0891 

.0685 

.0603 

.0560 

.0539 

.0533 

.0537 

.0550 

6 

.9651 

.7011 

.2839 

.1722 

.1459 

.1327 

.1233 

.1157 

.1093 

.1037 

.0987 

7 

.9646 

.7099 

.2448 

.0954 

.0623 

.0484 

.0394 

.0329 

.0278 

.0237 

.0203 

8 

.9752 

.8121 

.3666 

.1191 

.0604 

.0416 

.0320 

.0260 

.0220 

.0194 

.0178 

9 

.9705 

.8183 

.4487 

.2143 

.1427 

.1156 

.0999 

.0884 

.0793 

.0718 

.0654 

10 

.9775 

.8722 

.5502 

.2312 

.1085 

.0677 

.0488 

.0370 

.0287 

.0224 

.0177 

11 

.9691 

.8539 

.5671 

.2970 

.1771 

.1299 

.1061 

.0907 

.0793 

.0702 

.0628 

12 

.9792 

.9085 

.6966 

.3901 

.1955 

.1142 

.0786 

.0590 

.0460 

.0365 

.0293 

13 

.9729 

.8991 

.7126 

.4566 

.2738 

.1829 

.1383 

.1128 

.0956 

.0827 

.0725 

14 

.9816 

.9356 

.8095 

.5843 

.3608 

.2243 

.1556 

.1191 

.0968 

.0811 

.0692 

15 

.9838 

.9493 

.8583 

.6791 

.4535 

.2761 

.1733 

.1186 

.0879 

.0687 

.0555 

16 

.9851 

.9598 

.8993 

.7769 

.5904 

.3917 

.2401 

.1464 

.0930 

.0623 

.0439 

4.2  Dosimeter 

4.2.1  Motivation 

This  section  describes  the  modeling  of  the  CRRES  dosimeters  for  the  following  angular  distributions; 
the  isotropic  distribution,  the  mirror  plane  distribution  and  the  sin^  o  pitch  angle  distribution.  The 
mirror  plane  distribution  is  the  limiting  case  in  which  charged  particles  (in  this  case,  protons)  trapped 
in  the  earth’s  magnetic  field  are  at  their  "mirror  points”.  The  latter  are  the  points  at  which  the  pitch 
angle  (the  angle  between  the  magnetic  field  vector  and  the  proton’s  velocity  vector)  is  90  degrees.  The 
mirror  points  represent  the  northern  or  southern  extremes  of  their  motion  relative  to  the  geomagnetic 
field.  The  angular  distribution  of  mirroring  particles  is  uniform  in  the  mirror  plane,  and  is  sometimes 
referred  to  as  a  "delta  function"  in  the  mirror  plane.  The  software  models  for  these  angular  distributions 
will  be  used  as  an  aid  to  analyze  the  data  produced  by  the  CRRES  Dosimeters,  which,  in  turn,  will  be 
used  to  model  the  proton  distribution  in  the  radiation  belts. 


The  design  of  the  CRRES  Dosimeters  is  based  on  the  Space  Radiation  Dosimeters  successfully  used  on 
the  DMSP  7  satellite  IGussenhoven,  et.  al.,  1986|.  Modifications  of  the  original  DMSP  design  were 
made  because  of  the  different  operational  environment  and  the  orbit  of  the  CRRES  satellite.  A 
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description  of  the  CRRES  Dosimeters  may  be  found  in  Morel,  et.  al.  [1989].  An  overall  description 
of  the  CRRES  Satellite  and  Experiments  can  be  found  in  Gussenkoven,  et.  al.  [1987].  Each  of  the 
CRRES  Dosimeter  detectors  consists  of  a  hemispherical  aluminum  housing,  which  is  used  to  exclude 
charged  particles  with  energies  below  a  certain  threshold  limit.  At  the  center  of  the  hemispheric  shell 
is  a  silicon  disk-shaped  particle  detector  which  measures  the  energy  deposited  in  its  active  region  by  high 
energy  ions,  and  reports  the  counts  and  dose  as  follows:  there  are  three  ranges  (channels)  to  which  the 
instruments  respond,  viz,  LOLE'I'  (energy  depositions  of  0.5  -  1  MeV),  HILET  (energy  depositions  of 
1  -  10  MeV),  and  star  events  (energy  depositions  greater  than  20  MeV).  In  this  report,  we  will  only  be 
concerned  with  the  LOLET  and  HILET  channels,  and  with  protons  in  the  mirror  plane  particle 
distribution  limit. 

4.2.2  Methodology 

The  approach  used  to  model  the  CRRES  Dosimeters,  in  this  study,  is  to  compute  the  path  length 
distribution,  and  use  it  to  compute  the  counts  and  dose  for  a  variety  of  mirror  plane  proton  flux 
distributions.  Ilie  path  length  distribution  represents  a  "histogram"  of  the  possible  path  lengths  which 
a  charged  particle  (proton),  from  the  mirror  plane  distribution,  can  take  through  the  detector,  and  is  a 
function  of  the  angle  between  the  magnetic  field  vector  and  the  normal  to  the  top  surface  of  the  detector. 
If  normalized  to  1,  the  path  length  distribution  represents  the  probability  that  a  particle  passing  straight 
through  the  detector  will  have  a  given  path  length. 

The  path  length  distribution  depends  only  on  the  detector  geometry  and  on  the  angular  particle 
distribution.  Using  code  based  upon  the  Janni  Energy  -  Range  tables,  the  path  length  distribution  is 
used  to  compute  the  counts  and  dose  which  would  be  observed  for  a  given  proton  energy  spectrum. 


% 


4.2.2.1 


Physical  Assumptions 

The  CRRES  Dosimeters  are  a  set  of  four  disk  shaped  siiicort  particle  detectors,  each  of  which  is  mounted 
at  the  center  of  a  hemispheric  aluminum  shell.  For  simplicity  in  computation,  only  the  active  portion 
of  the  detectors  is  modeled.  The  hemispherical  shell  is  used  to  set  the  lower  energy  limits  of  charged 
particles  which  will  arrive  at  the  detector.  'Ihe  physical  and  geometrical  parameters  of  the  detectors  are 
presented  in  Table  6. 


Table  6.  Physical  and  Geometric  Parameters  of  the  CRRES  Dosimeter  Detectors 

Detector  # 

Detector  Area 
[cm^l 

Detector 

Thickness 

(microns] 

Thickness  of  AJ 
Hemisphere 
[gm/cm^j 

Minimum  Proton 
Energy  Required  to 
Penetrate  Aluminum 
Hemisphere  [MeV] 

1 

0.00815 

403 

0.55 

20 

2 

0.051 

434 

1..55 

35 

3 

0.051 

399 

3.05 

51 

4 

1.000 

406 

5.91 

75 

A  summary  of  the  notation  used  in  this  report  is  provided  in  the  following  table: 


Table  7.  Notation 

Symbol 

Description 

Units 

D 

Thickness  of  Detector 

microns 

R 

Radius  of  Detector 

cm 

X 

Angle  between  magnetic  field  lines  and  normal  to  detector  surface 

Y 

Angle  in  mirror  plane  (measured  from  the  projection  of  the  normal  to 
detector  surface  onto  the  mirror  plane) 

B 

Magnetic  Field  Vector 

n 

Unit  vector  normal  to  detector  top  surface 

The  detector  geometry  for  the  mirror  plane  case  is  illustrated  in  Figure  51. 
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O«t«ctori  Racfcj*  »  R,  Tbickn»*«  *  D  intersection  o-f  «  mirror  p)aoe  with  disk 


Figure  51.  Detector  Geometry 


The  mirror  plane  geometry  is  defined  by  X,  the  angle  between  the  magnetic  field  vector  B  and  the 
normal  to  the  detector  surface  n.  Each  plane  perpendicular  to  B  is  a  mirror  plane.  For  each  point, 
either  on  the  detector  top  surface,  or  within  the  volume  of  the  detector,  there  is  a  single  mirror  plane 
upon  which  the  point  lies.  In  the  coordinate  system  used  here,  the  vector  n  points  towards  the  positive 
z-axis,  and  the  x-axis  lies  in  the  plane  defined  by  B  and  n.  ITie  y-axis  is  perpendicular  to  both  B  and  n. 
The  angle  y  'S  defined  as  the  angle  in  the  mirror  plane  between  the  velocity  vector  and  the  projection 
of  the  vector  -n  onto  the  mirror  plane. 


4.2.22  Path  Length  Distribution  Computations 

The  problem  is,  therefore,  to  determine  the  path  length  distribution,  as  a  weighted  sum,  of  the 
contributions  from  different  directions  as  determined  by  the  geometry  of  the  detector  and  of  the  particle 
angular  distribution.  Monte  Carlo  methods  were  used  to  compute  the  path  length  distributions.  A  full 
Monte  Carlo  computation  was  numerically  impractical,  because  it  would  require  the  selection  of  large 
numbers  of  points  lying  outside  the  detector,  and  would  require  a  determination  for  each  trial  whether 
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the  ray  selected  would  actually  penetrate  the  detector.  It  is,  therefore,  desirable  to  compute  the 
contributions  due  to  the  following  possible  cases: 


(1)  A  proton  enters  the  top  of  the  detector,  and  exits  the  bottom  or  the  side  of  the  detector. 

(2)  A  proton  enters  the  side  of  the  detector,  and  exits  through  either  the  side  or  the  bottom. 

ITiese  Monte  Carlo  computations  will  be  briefly  described  below  for  the  Mirror  Plane  case. 

Case  1. 

The  Monte  Carlo  program  picks  random  points  P  on  the  top  surface  of  the  detector.  In  the  mirror 
plane  passing  through  P,  a  value  of  y  was  selected  randomly  according  to  the  proper  probability 
distribution.  A  trial  value  of  q  was  computed  using  the  equation  q  =  D/sin  X  cos  y-  The  vector  q  u 
from  P,  in  the  direction  defined  by  X  and  y,  was  computed.  If  the  resulting  point  O  =  (x',  y'.  z') 
satisfied  the  condition  x'^  +  y'^  <  n,  then  a  counter,  corresponding  to  an  interval  (bin)  for  the  value  of 
q  (top-to-bottom),  is  incremented.  Otherwise,  the  quadratic  equation  in  q. 


sin^  X  cos^  y)  +  2  {x  cos  y  cos  X  +  y  sin  y)  j  +  {x" 


-  R^)  =  0 


which  implements  the  conditions  x'^  +  y'^  =  0  <  z'  <  D,  is  solved,  and  the  counter  (top-to-side 

case),  corresponding  to  q,  is  incremented. 

Case  2. 

For  the  Monte  Carlo  computation,  points  are  randomly  selected  in  the  y-z  plane,  -R  <  y  <  R,  0  <  z 
<  D;  x  was  computed  from  the  equation  x  =  -  V(R^  -  y^),  and  y  was  computed  from  a  uniform 
distribution  on  the  interval  -n/2  <  y  <  Jr/2.  First,  a  trial  value  of  q  was  computed  using  the  equation  q 
=  z  /  sin  X  cos  y.  If  x'^  +  y'^  <  R^  then  the  side-bottom  counter  was  incremented  for  the  interval 
containing  q.  If  the  trial  fails,  then  q  is  computed  using  the  equation, 

q  =  -2(x  cos  y  cos  X  y  sin  y) 

1  -  sin^  X  cos^  y 


and  the  side-side  counter  was  incremented,  provided  that  z"  >  0  and  q  was  in  the  proper  range. 

To  combine  the  Case  1  and  Case  2  Monte  Carlo  results  into  a  single  path  length  distribution,  it  is 
necessary  to  determine  the  relative  flux  contributions  entering  from  the  top  and  from  the  side.  As  seen 
from  the  side,  the  area  of  the  detector  is  2  D  R.  If  we  multiply  the  area  by  cos  X.  we  obtain  the  projected 
area  as  seen  from  protons  in  the  mirroi  plane  incident  on  the  detector.  We  wish  to  normalize  the 
computation  with  respect  to  the  top  surface  area  of  the  detector,  which  is  given  by  nR^.  The  Case  2 
contribution  is  obtained  from  the  normalized  Case  2  distribution  by  multiplying  the  latter  by  2  D  cos 
X/  n  R.  Similarly,  the  projected  area  of  the  top  surfaa  as  seen  from  protons  entering  the  detector  from 
the  top  must  be  multiplied  sin  X.  It  is  also  necessary  to  take  into  account,  that  for  any  value  of  X, 
there  is  always  a  contribution  to  the  side-side  component  corresponding  to  y  =  ±rt/2,  while  for  X  «  0, 
the  side  to  side  contribution  is  present  for  each  value  of  y.  Therefore,  the  total  contribution  is  given  by: 
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f(>.)  -2D  (1/ji  +  cos  X  (1  -  l/ii))/R  S(X)  +  sin  X  T(X) 


where  S  and  T  are  the  contributions  from  the  side  and  from  the  top,  respectively.  The  function  f(X)  will 
then  represent  the  path  length  distribution  (relative  to  unit  area). 

4.2.3  Results 

The  actual  mirror  plane  path  length  distributions  were  computed  for  values  of  X  from  0  to  90  degrees 
at  one  degree  intervals.  Figures  52  and  53  represent  the  mirror  plane  path  length  distributions  for 
detectors  1  and  4  for  X  =  60°.  For  comparison,  the  infinite  slab  mirror  plane  approximation  path 
length  distribution  is  also  provided  on  the  same  graphs  using  the  same  normalization  (area  under  the 
curve).  Note  the  close  agreement  for  Figure  53;  for  Figure  52  the  differences  are  due  to  edge  effects. 
The  infinite  slab  path  length  distribution  is  zero  until  the  minimum  possible  path  length  which  occurs 
at  D/sm  X. 
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Figure  52.  Path  Length  Distribution  for  Detector  1  for  X  =  60  degrees  compared  with  Infinite 
Slab  Path  Length  Distribution  (smooth  curve). 
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Figure  53.  Path  Length  Distribution  for  Detector  4,  A.  =  60  degrees. 


To  compute  the  flux  and  dose  response  function,  for  each  incident  energy  value,  the  energy  loss  is 
computed  for  each  path  length  value,  and  the  corresponding  computed  count  and  dose  contribution  is 
weighted  by  the  normalized  path  length  distribution  value. 
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